[petsc-users] petsc vector management for monte carlo

[Matthew Knepley]

On Sat, Nov 30, 2013 at 2:56 PM, Gideon Simpson <gideon.simpson at gmail.com>wrote:

> Thanks Jed, see a few follow up comments below,
> On Nov 30, 2013, at 3:34 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
>
> > Gideon Simpson <gideon.simpson at gmail.com> writes:
> >
> >> I'm migrating some monte carlo code over to petsc, where the
> realizations will be done serially, but distributed across the cluster.
>  I'm inclined to involve petsc in this for its data management, and had a
> couple of questions.
> >>
> >> If I plant to store the (scalar) results of the monte carlo simulations
> in a Petsc vector v, what is the recommended method for putting the results
> into this vector?  Currently, I have implemented:
> >>
> >> VecGetOwnershipRange(sim_data, &local_start, &local_end);
> >>
> >> VecGetArray(sim_data,& sim_data_local);
> >> k=0;
> >>
> >> for(i=local_start;i<local_end;i++){
> >>
> >>    sim_result = mc_sim_cod(mc_args);
> >>
> >>    sim_data_local[k] = sim_result;
> >>    k++;
> >>
> >>  }
> >> VecRestoreArray(sim_data, & sim_data_local);
> >>
> >> Is this a *good* solution, or should I be using SetValue routines?
> >
> > The above is fine.
> >
> >> Also, I have MC simulations which, instead of returning a single
> >> scalar, return some n-dimensional array of real values.  In that case,
> >> I was thinking to make sim_data a petsc matrix.  If each result was
> >> n-dimensional, how would I initialize the matrix to hold M
> >> realizations of this n dimensional matrix?  ensuring that each set of
> >> n was contiguous on the local machine?
> >
> > Mat is really for linear operators, not a container for data that you
> > want to interpret as 2-dimensional.  DMDA can do the latter, especially
> > if you want to distribute.
> >
> >> Alternatively, this could be an n*M vector, but again, I'd need to
> >> ensure each chunk of n was located on the same machine?
> >
> > Just set the local sizes of the Vec such that the entire subset is
> > local.  See the middle argument to VecSetSizes().
>
> If each realization generated a size n output, is there something slicker
> than this:
>
> VecGetOwnershipRange(sim_data, &local_start, &local_end);
>
> VecGetArray(sim_data,& sim_data_local);
> k=0;
>
> m_local = (local_end - local_start-1)/n;
>
> for(i=0;i< m_local;i++){
>
>    mc_sim_cod(mc_args, &sim_result);
>
>    for(k = 0;k<n;k++){
>
>       sim_data_local[i * n + k] = sim_result[k];
>    }
>
>  }
> VecRestoreArray(sim_data, & sim_data_local);
>
> I see  expressions like sim_data_local[i * n + k] to handle the indexing
> as  bit awkward.
>

I guess it depends on what you want to do. You could have n independent
serial vectors.
If you want to maintain the global vector, you can always define the
pointer at the top of
the loop.

  Thanks,

     Matt


> >
> >> I presume here, I would use SetValue routines to insert each chunk of
> >> data.
> >
> > You could, but writing directly into the array is also fine.
> >
> >> Also, this may have to do with it being the newest matlab, 2013b, but
> when I try to load petsc vectors i receive the following error:
> >>
> >> Error using load
> >> Number of columns on line 2 of ASCII file sim_data.out
> >> must be the same as previous lines.
> >>
> >> where I constructed the sim_data.out file using the lines:
> >>
> >> PetscViewerASCIIOpen(PETSC_COMM_WORLD, "sim_data.out", &output_viewer);
> >> PetscViewerSetFormat(output_viewer, PETSC_VIEWER_ASCII_MATLAB);
> >> VecView(sim_data,output_viewer);
> >> PetscViewerDestroy(&output_viewer);
> >>
> >> Looking at the ascii file, it looks fine.
> >
> > Did MATLAB change their ASCII format?
>
> Here's a quick test on 2013b:
> v = linspace(0,1,5);
> save('test.out', 'v', '-ascii');
>
> test.out looks like:
>    0.0000000e+00   2.5000000e-01   5.0000000e-01   7.5000000e-01
> 1.0000000e+00
>
>
> >
> > Anyway, it is much better to write binary output and read with
> > PetscBinaryRead (bin/matlab/PetscBinaryRead.m).
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131130/0f0455eb/attachment.html>