[petsc-users] petsc vector management for monte carlo

[Gideon Simpson]

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 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).