Matthew Knepley wrote: > On Mon, Jul 5, 2010 at 3:22 PM, Thomas Witkowski > <thomas.witkowski at tu-dresden.de > <mailto:thomas.witkowski at tu-dresden.de>> wrote: > > Hi, > > I've some trouble with matrix values that are set by MatSetValues, > but are missing in the final matrix. I can reduce the problem to a > 75x75 matrix that is created on four processors. I create it quite > simple with: > > MatCreateMPIAIJ(PETSC_COMM_WORLD, 12, 12, 75, 75, petscMatrix) on p0 > MatCreateMPIAIJ(PETSC_COMM_WORLD, 18, 18, 75, 75, petscMatrix) on p1 > MatCreateMPIAIJ(PETSC_COMM_WORLD, 18, 18, 75, 75, petscMatrix) on p2 > MatCreateMPIAIJ(PETSC_COMM_WORLD, 27, 27, 75, 75, petscMatrix) on p3 > > On all processors, the values are set with the following command: > > MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), > &(values[0]), ADD_VALUES); > > where rowIndex is an integer, cols is of type std::vector<int> and > values if of type std::vector<double>. Before > > > I do not think this is always kosher. std::vector<> may store values > however it wants (perhaps in a > contiguous array). std::valarray<> I believe mandates a contiguous > array, however even it does not > give you access to the explicit pointer (you have to use &v[0] > nonsense). Reasons to hate the STL. This is guaranteed by the C++ standard, so this shouldn't be the problem:
A vector's storage is guaranteed to be contiguous. From ISO 14882, 2nd ed., 23.2.4 [lib.vector]: "The elements of a vector are stored contiguously, meaning that if v is a vector<T, Allocator> where T is some type other than bool, then it obeys the identity &v[n] == &v[0] + n for all 0 <= n < v.size()." > > > MatSetValues is called, I run over the arrays and print all the > entries that are added to the matrix. Finally, > > MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY); > MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY); > > are called. To check the matrix, I use the option > -mat_view_matlab. Okay, now some of the entries are missing and > I've absolute no idea what I did wrong (i.e. entry row 48-col 69, > which is only once by rank 3). What is the best way to debug the > problem? I use petsc 3.0.0p11, which is compiled in debug mode. My > code is also compiled in debug mode and valgrind does not report > any errors. Thanks for any suggestions. > > > MatAssemble and MatView after every insertion. Good idea, will try it. Thomas > > Matt > > > > Thomas > > > > > -- > 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
