Send the code to petsc-maint at mcs.anl.gov likely you are not calling MatAssemblyBegin/End or something similar in the correct place.
Barry On Dec 17, 2010, at 9:23 AM, Mohamad M. Nasr-Azadani 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! > 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 > > >
