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