[petsc-users] How to debug problem with MatSetValues

[Thomas Witkowski]

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