[petsc-users] Matrix setup and Distributed arrays,

[Matthew Knepley]

On Fri, Dec 17, 2010 at 7:23 AM, Mohamad M. Nasr-Azadani
<mmnasr at gmail.com>wrote:

> 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!
>

I do not know what you mean by the "right result". If the only difference is
a 0 in the matrix, you will
get the same result for MatMult(). Did you check this?

  Thanks,

     Matt


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


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20101217/ef27cf9e/attachment.htm>