Program Listing for File bigfeta_dist.h

↰ Return to documentation for file (bigfeta/distributed/src/bigfeta_dist.h)

#include <petscsys.h>
#include <petscksp.h>
#include <stdio.h>
#include <stdlib.h>
#include <petscviewerhdf5.h>
#include <libgen.h>
#include <hdf5.h>

PetscErrorCode
CountFiles (MPI_Comm COMM, char indexname[], int *nfiles)
{
  PetscErrorCode ierr;
  PetscMPIInt rank;
  ierr = MPI_Comm_rank (COMM, &rank);
  CHKERRQ (ierr);

  if (rank == 0)
    {
      hid_t file, space, dset;
      hsize_t dims[1];

      file = H5Fopen (indexname, H5F_ACC_RDONLY, H5P_DEFAULT);
      dset = H5Dopen (file, "datafile_names", H5P_DEFAULT);
      space = H5Dget_space (dset);
      H5Sget_simple_extent_dims (space, dims, NULL);
      *nfiles = dims[0];
      H5Dclose (dset);
      H5Sclose (space);
      H5Fclose (file);
    }
  MPI_Bcast (nfiles, 1, MPI_INT, 0, COMM);
  return ierr;
}

PetscErrorCode
ReadMetadata (MPI_Comm COMM, char indexname[], int nfiles, char *csrnames[],
          PetscInt ** metadata)
{
  PetscErrorCode ierr;
  PetscMPIInt rank;
  int i;
  ierr = MPI_Comm_rank (COMM, &rank);
  CHKERRQ (ierr);

  if (rank == 0)
    {
      char **rdata;
      hid_t file, filetype, memtype, space, dset;
      hsize_t dims[1];
      char *dir, *tmp;

      tmp = strdup (indexname);
      dir = strdup (dirname (tmp));
      file = H5Fopen (indexname, H5F_ACC_RDONLY, H5P_DEFAULT);
      dset = H5Dopen (file, "datafile_names", H5P_DEFAULT);
      filetype = H5Dget_type (dset);
      space = H5Dget_space (dset);
      H5Sget_simple_extent_dims (space, dims, NULL);
      rdata = (char **) malloc (dims[0] * sizeof (char *));
      memtype = H5Tcopy (H5T_C_S1);
      H5Tset_size (memtype, H5T_VARIABLE);
      H5Dread (dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, rdata);
      for (i = 0; i < nfiles; i++)
    {
      sprintf (csrnames[i], "%s/%s", dir, rdata[i]);
    }
      H5Dvlen_reclaim (memtype, space, H5P_DEFAULT, rdata);
      free (rdata);
      H5Dclose (dset);
      H5Sclose (space);

      PetscInt *row, *mxcol, *mncol, *nnz;
      row = (PetscInt *) malloc (nfiles * sizeof (PetscInt));
      mxcol = (PetscInt *) malloc (nfiles * sizeof (PetscInt));
      mncol = (PetscInt *) malloc (nfiles * sizeof (PetscInt));
      nnz = (PetscInt *) malloc (nfiles * sizeof (PetscInt));
      dset = H5Dopen (file, "datafile_nrows", H5P_DEFAULT);
      H5Dread (dset, H5T_NATIVE_LONG, H5S_ALL, H5S_ALL, H5P_DEFAULT, row);
      dset = H5Dopen (file, "datafile_mincol", H5P_DEFAULT);
      H5Dread (dset, H5T_NATIVE_LONG, H5S_ALL, H5S_ALL, H5P_DEFAULT, mncol);
      dset = H5Dopen (file, "datafile_maxcol", H5P_DEFAULT);
      H5Dread (dset, H5T_NATIVE_LONG, H5S_ALL, H5S_ALL, H5P_DEFAULT, mxcol);
      dset = H5Dopen (file, "datafile_nnz", H5P_DEFAULT);
      H5Dread (dset, H5T_NATIVE_LONG, H5S_ALL, H5S_ALL, H5P_DEFAULT, nnz);
      for (i = 0; i < nfiles; i++)
    {
      metadata[i][0] = row[i];
      metadata[i][1] = mncol[i];
      metadata[i][2] = mxcol[i];
      metadata[i][3] = nnz[i];
    }
      free (row);
      free (mncol);
      free (mxcol);
      free (nnz);
      H5Dclose (dset);
      H5Tclose (filetype);
      H5Tclose (memtype);
      H5Fclose (file);
      free (dir);
      free (tmp);
    }
  for (i = 0; i < nfiles; i++)
    {
      MPI_Bcast (csrnames[i], PETSC_MAX_PATH_LEN, MPI_CHAR, 0, COMM);
      MPI_Bcast (metadata[i], 4, MPIU_INT, 0, COMM);
    }
  return ierr;
}

void
CopyDataSetstoSolutionOut (MPI_Comm COMM, char indexname[], char outputname[])
{
  hid_t filein, fileout, filetype, memtype, space, dset, dsetout;
  hsize_t dims[1];
  int nds = 8;
  const char *copyids[8] = { "input_args", "resolved_tiles",
    "datafile_names", "datafile_maxcol", "datafile_mincol",
    "datafile_nnz", "datafile_nrows", "reg", "solve_list"
  };
  char **rdata;
  int i;

  filein = H5Fopen (indexname, H5F_ACC_RDONLY, H5P_DEFAULT);
  fileout = H5Fcreate (outputname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
  for (i = 0; i < nds; i++)
    {
      dset = H5Dopen (filein, copyids[i], H5P_DEFAULT);
      filetype = H5Dget_type (dset);
      space = H5Dget_space (dset);
      H5Sget_simple_extent_dims (space, dims, NULL);
      rdata = (char **) malloc (dims[0] * sizeof (char *));
      memtype = H5Dget_type (dset);
      if (i < 3)
    {
      memtype = H5Tcopy (H5T_C_S1);
      H5Tset_size (memtype, H5T_VARIABLE);
    }
      H5Dread (dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, rdata);
      dsetout =
    H5Dcreate (fileout, copyids[i], filetype, space, H5P_DEFAULT,
           H5P_DEFAULT, H5P_DEFAULT);
      H5Dwrite (dsetout, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, rdata);
      H5Dclose (dset);
      H5Dclose (dsetout);
      H5Tclose (filetype);
      H5Tclose (memtype);
      H5Dvlen_reclaim (memtype, space, H5P_DEFAULT, rdata);
      free (rdata);
    }
  H5Fclose (filein);
  H5Fclose (fileout);
}

PetscErrorCode
SetFiles (MPI_Comm COMM, int nfiles, PetscInt * firstfile,
      PetscInt * lastfile)
{
  PetscErrorCode ierr;
  PetscMPIInt rank, size;
  float avg;
  int i, rem, last;

  ierr = MPI_Comm_rank (COMM, &rank);
  CHKERRQ (ierr);
  ierr = MPI_Comm_size (COMM, &size);
  CHKERRQ (ierr);
  if (nfiles <= size)
    {
      for (i = 0; i < size; i++)
    {
      if (rank == i)
        {
          *firstfile = i;
          if (i > nfiles - 1)
        {
          *lastfile = i - 1;
        }
          else
        {
          *lastfile = i;
        }
        }
    }
    }
  else
    {
      avg = nfiles / (float) size;
      rem = nfiles % size;
      last = 0;
      for (i = 0; i < size; i++)
    {
      if (rank == i)
        {
          *firstfile = last;
        }
      if (rem-- > 0)
        {
          last = last + ceil (avg);
        }
      else
        {
          last = last + floor (avg);
        }
      if (rank == i)
        {
          *lastfile = last - 1;
        }
    }
    }
  return ierr;
}

PetscErrorCode
ReadVec (MPI_Comm COMM, PetscViewer viewer, char *varname, Vec * newvec,
     PetscInt * n)
{
  PetscErrorCode ierr;
  char name[256];
  sprintf (name, "%s", varname);
  ierr = VecCreate (COMM, newvec);
  CHKERRQ (ierr);
  ierr = PetscObjectSetName ((PetscObject) * newvec, name);
  CHKERRQ (ierr);
  ierr = VecLoad (*newvec, viewer);
  CHKERRQ (ierr);
  ierr = VecGetSize (*newvec, n);
  CHKERRQ (ierr);
  return ierr;
}

PetscErrorCode
ReadVecWithSizes (MPI_Comm COMM, PetscViewer viewer, char *varname,
          Vec * newvec, PetscInt * n, PetscInt nlocal,
          PetscInt nglobal, PetscBool trunc)
{
  PetscErrorCode ierr;
  Vec tmpvec;
  PetscInt i, ntmp, low, high, *indx;
  PetscScalar *vtmp;
  char name[256];

  sprintf (name, "%s", varname);
  ierr = VecCreate (COMM, newvec);
  CHKERRQ (ierr);
  ierr = VecSetSizes (*newvec, nlocal, nglobal);
  CHKERRQ (ierr);
  ierr = VecSetType (*newvec, VECMPI);
  CHKERRQ (ierr);
  if (trunc)
    {
      ierr = VecCreate (MPI_COMM_SELF, &tmpvec);
      CHKERRQ (ierr);
      ierr = VecSetType (tmpvec, VECSEQ);
      CHKERRQ (ierr);
      ierr = PetscObjectSetName ((PetscObject) tmpvec, name);
      CHKERRQ (ierr);
      ierr = VecLoad (tmpvec, viewer);
      CHKERRQ (ierr);
      ierr = VecGetOwnershipRange (*newvec, &low, &high);
      CHKERRQ (ierr);
      indx = (PetscInt *) malloc ((high - low) * sizeof (PetscInt));
      vtmp = (PetscScalar *) malloc ((high - low) * sizeof (PetscScalar));
      for (i = low; i < high; i++)
    {
      indx[i - low] = i;
    }
      ierr = VecGetValues (tmpvec, high - low, indx, vtmp);
      CHKERRQ (ierr);
      ierr = VecSetValues (*newvec, high - low, indx, vtmp, INSERT_VALUES);
      CHKERRQ (ierr);
      free (indx);
      free (vtmp);
      ierr = VecAssemblyBegin (*newvec);
      CHKERRQ (ierr);
      ierr = VecAssemblyEnd (*newvec);
      CHKERRQ (ierr);
      ierr = VecDestroy (&tmpvec);
      CHKERRQ (ierr);
    }
  else
    {
      ierr = PetscObjectSetName ((PetscObject) * newvec, name);
      CHKERRQ (ierr);
      ierr = VecLoad (*newvec, viewer);
      CHKERRQ (ierr);
    }
  ierr = VecGetSize (*newvec, &ntmp);
  CHKERRQ (ierr);
  return ierr;
}

PetscErrorCode
ReadIndexSet (MPI_Comm COMM, PetscViewer viewer, char *varname, IS * newIS,
          PetscInt * n)
{
  PetscErrorCode ierr;
  char name[256];
  sprintf (name, "%s", varname);
  ierr = ISCreate (COMM, newIS);
  CHKERRQ (ierr);
  ierr = PetscObjectSetName ((PetscObject) * newIS, name);
  CHKERRQ (ierr);
  ierr = ISLoad (*newIS, viewer);
  CHKERRQ (ierr);
  ierr = ISGetSize (*newIS, n);
  CHKERRQ (ierr);
  return ierr;
}

PetscErrorCode
ShowMatInfo (MPI_Comm COMM, Mat * m, const char *mesg)
{
  PetscErrorCode ierr;
  MatInfo info;
  PetscBool isassembled;
  PetscInt rowcheck, colcheck;
  PetscMPIInt rank;

  ierr = MPI_Comm_rank (COMM, &rank);
  CHKERRQ (ierr);
  ierr = MatGetInfo (*m, MAT_GLOBAL_SUM, &info);
  CHKERRQ (ierr);
  if (rank == 0)
    {
      printf ("%s info from rank %d:\n", mesg, rank);
      ierr = MatAssembled (*m, &isassembled);
      CHKERRQ (ierr);
      printf (" is assembled: %d\n", isassembled);
      ierr = MatGetSize (*m, &rowcheck, &colcheck);
      CHKERRQ (ierr);
      printf (" global size %ld x %ld\n", rowcheck, colcheck);
      //ierr = MatGetInfo(*m,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
      printf (" block_size: %f\n", info.block_size);
      printf (" nz_allocated: %f\n", info.nz_allocated);
      printf (" nz_used: %f\n", info.nz_used);
      printf (" nz_unneeded: %f\n", info.nz_unneeded);
      printf (" memory: %f\n", info.memory);
      printf (" assemblies: %f\n", info.assemblies);
      printf (" mallocs: %f\n", info.mallocs);
      printf (" fill_ratio_given: %f\n", info.fill_ratio_given);
      printf (" fill_ratio_needed: %f\n", info.fill_ratio_needed);
    }
  ierr = MatGetInfo (*m, MAT_LOCAL, &info);
  CHKERRQ (ierr);
  if (rank == 0)
    {
      printf ("%s local info from rank %d:\n", mesg, rank);
      ierr = MatAssembled (*m, &isassembled);
      CHKERRQ (ierr);
      printf (" is assembled: %d\n", isassembled);
      ierr = MatGetSize (*m, &rowcheck, &colcheck);
      CHKERRQ (ierr);
      printf (" global size %ld x %ld\n", rowcheck, colcheck);
      //ierr = MatGetInfo(*m,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr);
      printf (" block_size: %f\n", info.block_size);
      printf (" nz_allocated: %f\n", info.nz_allocated);
      printf (" nz_used: %f\n", info.nz_used);
      printf (" nz_unneeded: %f\n", info.nz_unneeded);
      printf (" memory: %f\n", info.memory);
      printf (" assemblies: %f\n", info.assemblies);
      printf (" mallocs: %f\n", info.mallocs);
      printf (" fill_ratio_given: %f\n", info.fill_ratio_given);
      printf (" fill_ratio_needed: %f\n", info.fill_ratio_needed);
    }
  return ierr;
}

void
GetGlobalLocalCounts (int nfiles, PetscInt ** metadata, int local_firstfile,
              int local_lastfile, PetscInt * global_nrow,
              PetscInt * global_ncol, PetscInt * global_nnz,
              PetscInt * local_nrow, PetscInt * local_nnz,
              PetscInt * local_row0)
{
  PetscInt cmin, cmax;
  int i;

  *global_nrow = 0;
  *global_ncol = 0;
  *global_nnz = 0;
  *local_nrow = 0;
  *local_nnz = 0;
  *local_row0 = 0;
  cmin = 0;
  cmax = 0;
  for (i = 0; i < nfiles; i++)
    {
      *global_nrow += metadata[i][0];
      *global_nnz += metadata[i][3];
      if (i < local_firstfile)
    {
      *local_row0 += metadata[i][0];
    }
      if (i == 0)
    {
      cmin = metadata[i][1];
      cmax = metadata[i][2];
    }
      else
    {
      if (metadata[i][1] < cmin)
        {
          cmin = metadata[i][1];
        }
      if (metadata[i][2] > cmax)
        {
          cmax = metadata[i][2];
        }
    }
    }
  *global_ncol = cmax - cmin + 1;
  for (i = local_firstfile; i <= local_lastfile; i++)
    {
      *local_nnz += metadata[i][3];
      *local_nrow += metadata[i][0];
    }
  return;
}

PetscErrorCode
ReadLocalCSR (MPI_Comm COMM, char *csrnames[], int local_firstfile,
          int local_lastfile, int nsolve, PetscInt * local_indptr,
          PetscInt * local_jcol, PetscScalar * local_data,
          PetscScalar * local_weights, PetscScalar ** local_rhs)
{

  PetscViewer viewer;       //viewer object for reading files
  IS indices, indptr;
  Vec data, weights, rhs;
  PetscErrorCode ierr;
  int i, j;
  char tmp[200];

  //create the distributed matrix A
  PetscInt vcnt, niptr;
  PetscInt roff = 0, roff2 = 0, zoff = 0;
  const PetscInt *iptr, *jcol;
  PetscInt poff = 0;
  PetscScalar *a, *w;
  //these will concat the multiple files per processor
  for (i = local_firstfile; i <= local_lastfile; i++)
    {
      //open the file
      ierr = PetscViewerHDF5Open (COMM, csrnames[i], FILE_MODE_READ, &viewer);
      CHKERRQ (ierr);
      //indptr
      ierr = ReadIndexSet (COMM, viewer, (char *) "indptr", &indptr, &niptr);
      CHKERRQ (ierr);
      ISGetIndices (indptr, &iptr);
      for (j = 1; j < niptr; j++)
    {
      local_indptr[j + roff] = iptr[j] + poff;
    }
      ISRestoreIndices (indptr, &iptr);
      poff = local_indptr[niptr - 1 + roff];
      roff += niptr - 1;

      //indices
      ierr = ReadIndexSet (COMM, viewer, (char *) "indices", &indices, &vcnt);
      CHKERRQ (ierr);
      ISGetIndices (indices, &jcol);
      memcpy (&local_jcol[zoff], jcol, vcnt * sizeof (PetscInt));
      ISRestoreIndices (indices, &jcol);

      //data
      ierr = ReadVec (COMM, viewer, (char *) "data", &data, &vcnt);
      CHKERRQ (ierr);
      VecGetArray (data, &a);
      memcpy (&local_data[zoff], a, vcnt * sizeof (PetscScalar));
      VecRestoreArray (data, &a);
      zoff += vcnt;

      //weights
      ierr = ReadVec (COMM, viewer, (char *) "weights", &weights, &vcnt);
      CHKERRQ (ierr);
      VecGetArray (weights, &w);
      memcpy (&local_weights[roff2], w, vcnt * sizeof (PetscScalar));
      VecRestoreArray (weights, &w);

      //rhs
      for (j = 0; j < nsolve; j++)
    {
      sprintf (tmp, "rhs_%d", j);
      ierr = ReadVec (COMM, viewer, tmp, &rhs, &vcnt);
      CHKERRQ (ierr);
      VecGetArray (rhs, &w);
      memcpy (&local_rhs[j][roff2], w, vcnt * sizeof (PetscScalar));
      VecRestoreArray (rhs, &w);
    }
      roff2 += vcnt;

      ierr = PetscViewerDestroy (&viewer);
      CHKERRQ (ierr);
      ISDestroy (&indptr);
      VecDestroy (&data);
      VecDestroy (&weights);
      ISDestroy (&indices);
    }

  return ierr;
}

PetscErrorCode
CreateW (MPI_Comm COMM, PetscScalar * local_weights, PetscInt local_nrow,
     PetscInt local_row0, PetscInt global_nrow, Mat * W)
{
  PetscErrorCode ierr;
  Vec global_weights;
  PetscInt *indx;
  int i;

  ierr = VecCreate (COMM, &global_weights);
  CHKERRQ (ierr);
  ierr = VecSetSizes (global_weights, local_nrow, global_nrow);
  CHKERRQ (ierr);
  ierr = VecSetType (global_weights, VECMPI);
  CHKERRQ (ierr);
  indx = (PetscInt *) malloc (local_nrow * sizeof (PetscInt));
  for (i = 0; i < local_nrow; i++)
    {
      indx[i] = local_row0 + i;
    }
  ierr =
    VecSetValues (global_weights, local_nrow, indx, local_weights,
          INSERT_VALUES);
  CHKERRQ (ierr);
  free (indx);

  ierr = MatCreate (PETSC_COMM_WORLD, W);
  CHKERRQ (ierr);
  ierr = MatSetSizes (*W, local_nrow, local_nrow, global_nrow, global_nrow);
  CHKERRQ (ierr);
  ierr = MatSetType (*W, MATMPIAIJ);
  CHKERRQ (ierr);
  ierr = MatMPIAIJSetPreallocation (*W, 1, NULL, 0, NULL);
  CHKERRQ (ierr);
  ierr = MatDiagonalSet (*W, global_weights, INSERT_VALUES);
  CHKERRQ (ierr);

  ierr = VecDestroy (&global_weights);
  CHKERRQ (ierr);
  return ierr;
}

PetscErrorCode
CreateL (MPI_Comm COMM, char indexname[], PetscInt local_nrow,
     PetscInt global_nrow, PetscBool trunc, Mat * L)
{
  PetscErrorCode ierr;
  PetscViewer viewer;
  Vec global_reg;
  PetscMPIInt rank;
  PetscInt junk;

  ierr = MPI_Comm_rank (COMM, &rank);
  CHKERRQ (ierr);

  ierr = PetscViewerHDF5Open (COMM, indexname, FILE_MODE_READ, &viewer);
  CHKERRQ (ierr);
  ierr =
    ReadVecWithSizes (COMM, viewer, (char *) "reg", &global_reg, &junk,
              local_nrow, global_nrow, trunc);
  CHKERRQ (ierr);

  ierr = MatCreate (PETSC_COMM_WORLD, L);
  CHKERRQ (ierr);
  ierr = MatSetSizes (*L, local_nrow, local_nrow, global_nrow, global_nrow);
  CHKERRQ (ierr);
  ierr = MatSetType (*L, MATMPIAIJ);
  CHKERRQ (ierr);
  ierr = MatMPIAIJSetPreallocation (*L, 1, NULL, 0, NULL);
  CHKERRQ (ierr);
  ierr = MatDiagonalSet (*L, global_reg, INSERT_VALUES);
  CHKERRQ (ierr);
  MatAssemblyBegin (*L, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd (*L, MAT_FINAL_ASSEMBLY);
  ierr = PetscViewerDestroy (&viewer);
  CHKERRQ (ierr);
  VecDestroy (&global_reg);
  return ierr;
}

PetscErrorCode
CountSolves (MPI_Comm COMM, char indexname[], PetscInt * nsolve)
{
  PetscErrorCode ierr;
  PetscViewer viewer;
  PetscInt junk;
  IS test;

  *nsolve = 0;
  ierr = PetscViewerHDF5Open (COMM, indexname, FILE_MODE_READ, &viewer);
  CHKERRQ (ierr);
  ierr = ReadIndexSet (COMM, viewer, (char *) "solve_list", &test, &junk);
  CHKERRQ (ierr);
  *nsolve = junk;
  ierr = PetscViewerDestroy (&viewer);
  CHKERRQ (ierr);
  return ierr;
}

PetscErrorCode
Readx0 (MPI_Comm COMM, char indexname[], PetscInt local_nrow,
    PetscInt global_nrow, PetscInt nsolve, PetscBool trunc, Vec x0[])
{
  PetscErrorCode ierr;
  PetscViewer viewer;
  PetscInt junk;
  char tmp[200];
  int i;

  ierr = PetscViewerHDF5Open (COMM, indexname, FILE_MODE_READ, &viewer);
  CHKERRQ (ierr);
  for (i = 0; i < nsolve; i++)
    {
      sprintf (tmp, "x_%d", i);
      ierr =
    ReadVecWithSizes (COMM, viewer, tmp, &x0[i], &junk, local_nrow,
              global_nrow, trunc);
      CHKERRQ (ierr);
    }
  ierr = PetscViewerDestroy (&viewer);
  CHKERRQ (ierr);
  return ierr;
}