.. _program_listing_file_EMaligner_distributed_src_ema.h:

Program Listing for File ema.h
==============================

|exhale_lsh| :ref:`Return to documentation for file <file_EMaligner_distributed_src_ema.h>` (``EMaligner/distributed/src/ema.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #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;
   }