[petsc-users] How to debug problem with MatSetValues

[Thomas Witkowski]

Jed Brown wrote:
> On Mon, 05 Jul 2010 15:22:50 +0200, Thomas Witkowski <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 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).  
>>     
>
> Are you checking (49,70) in Matlab (PETSc uses 0-based indexing)?  When
> you say the value is missing, is it equal to zero, or not existent in
> the sparsity pattern?
>   
Thanks Jed, I`ve totally forgotten that Matlab indices are 1-based! That 
was the whole problem.

Thomas

> Can you confirm that MatSetValues is really being called with
> rowIndex=48 and cols containing 69 (e.g. by setting a breakpoint there)?
>
>
> Matt, std::vector is guaranteed to be contiguous in revisions of C++98
> and subsequent draft standards.  I think all implementations observe
> this.
>
>   http://stackoverflow.com/questions/247738/is-it-safe-to-assume-that-stl-vector-storage-is-always-contiguous
>
> Jed
>
>