On Sat, Nov 30, 2013 at 2:56 PM, Gideon Simpson <[email protected]>wrote:
> Thanks Jed, see a few follow up comments below, > On Nov 30, 2013, at 3:34 PM, Jed Brown <[email protected]> wrote: > > > Gideon Simpson <[email protected]> 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
