[petsc-users] Matrix setup and Distributed arrays,

[Mohamad M. Nasr-Azadani]

Hi guys,

I am trying to solve a simpel Poisson equation on regular grid (test case
for my code).
This is my problem. Sorry for long discussion. I want to make it clear what
I am doing.


1- Create distributed array (same number of processors and parallel layout):
DA_3D_STAR: DA_STAR_STENCIL, width=1
2- Create distributed array (same number of processors and parallel layout):
DA_3D_BOX  : DA_BOX_STENCIL, width=3
3- Create Matrix (A) using:     ierr = DAGetMatrix(DA_3D_STAR, MATMPIAIJ,
&A);
4- Setup Matrix A in the corresponding fashion

       4-1    For all the regular local nodes (based on DA_3D_STAR), simply
use MatSetValuesStencil() to insert maximum 7-nonzerons in each row.
       4-2    For some special nodes (still local though) use a 9-point
interpolation equation. Theses nonzeros might extend outside of 1-layer of
ghost nodes, also it might be in the direction not corresponding to the
7-point stencil. Hence, it would be possible that some rows in the matrix
owns nonzeros not within the local+ghost node regions. MatSetValuesStencil()
does not work.
       Here, the DA_3D_BOX comes useful.
       4-3     Return the global index of any node within the range of
DA_3D_STAR, using the command
                 ierr = DAGetGlobalIndices(DA_3D_BOX, &nt, &global_indices);

       4-4     Now, I can freely use MatSetValues() to insert the 9-point
interpolation nonzeros into the matrix (using the returned global indices).

I thought this should work since MatSetValues() does not care if it is
inserting in the local part or the global section of the matrix on the
neighboring processor.
However, for a very simple test case (solution of Poisson equation), this
does not work in parallel.
I narrowed it down and came to a very simple but annoying conclusion. I
don't think the test code that I have suffers from any bugs.
But this is what I got,
When a new nonzero is inserted in (even) one row of the matrix where the new
nonzero corresponds to a node from a neighboring processor and outside the
width=1, STAR_STENCIL layout, using MatSetValues(), I get wrong results.
Even if I insert a "zero" value, but at the given place, unfortunately I get
wrong result.

Note that, for a matrix of 10*10*10, I even tested that for one single row
of a matrix, and unbelievably, I get wrong results.
I printed the matrix for two cases where there is a new (nonzero) in the
row, and there is not!
Comparing the two matrices and using

>>diff Bad Good

This is the only existing difference:

150c150
< row 149: (100, 0)  (148, 0)  (149, 1)  (150, 0)  (156, -1)  (198, 0)
(212, 0)
---
> row 149: (100, 0)  (148, 0)  (149, 1)  (150, 0)  (156, -1)  (198, 0)


They are exactly the same, except that there is new zero at (212) which
should not alter the results!
But I don't get the right result for the first case!
This is driving me crazy! I don't have any clue why this happens.

I can send you the test code if you think it helps.
But to me, it sounds like when the matrix is created using DAGetMatrix(),
there might be some sort of restriction to adding new nonzeros to the
locations not defined based on the stencil and not within the range of the
DA.

In advance, thank you so much.
Mohamad
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101217/dcbc6ae9/attachment.htm>