On Fri, Dec 17, 2010 at 7:23 AM, Mohamad M. Nasr-Azadani <[email protected]>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