[petsc-users] Problem with MatZeroRowsColumnsIS()
Matthew Knepley
knepley at gmail.com
Wed May 14 07:24:49 CDT 2014
On Tue, May 13, 2014 at 9:00 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> Given that Jed wrote MatZeroRowsColumns_MPIAIJ() it is unlikely to be
> wrong.
>
> You wrote
>
> Calculating the norm2 of the residuals defined above in each case gives:
> MatZeroRowsIS() 1cpu: norm(res,2) = 0
> MatZeroRowsIS() 4cpu: norm(res,2) = 0
> MatZeroRowsColumnsIS() 1cpu: norm(res,2) = 1.6880e-10
> MatZeroRowsColumnsIS() 4cpu: norm(res,2) = 7.3786e+06
>
I think the issue is that the RHS changes between 1 and 4 processes. This
could be a bug, but I have
gone over the code and our test which look correct. I think it could also
be a misintepretation of how
it works since you are using a composite matrix. Does it work if you put
everything in a single matrix?
Thanks,
Matt
> why do you conclude this is wrong? MatZeroRowsColumnsIS() IS suppose to
> change the right hand side in a way different than MatZeroRowsIS().
>
> Explanation. For simplicity reorder the matrix rows/columns so that
> zeroed ones come last and the matrix is symmetric. Then you have
>
> ( A B ) (x_A) = (b_A)
> ( B D ) (x_B) (b_B)
>
> with MatZeroRows the new system is
>
> ( A B ) (x_A) = (b_A)
> ( 0 I ) (x_B) (x_B)
>
> it has the same solution as the original problem with the give x_B
>
> with MatZeroRowsColumns the new system is
>
> ( A 0 ) (x_A) = (b_A) - B*x_B
> ( 0 I ) (x_B) (x_B)
>
> note the right hand side needs to be changed so that the new problem has
> the same solution.
>
> Barry
>
>
> On May 9, 2014, at 9:50 AM, Püsök, Adina-Erika <puesoek at uni-mainz.de>
> wrote:
>
> > Yes, I tested the implementation with both MatZeroRowsIS() and
> MatZeroRowsColumnsIS(). But first, I will be more explicit about the
> problem I was set to solve:
> >
> > We have a Dirichlet block of size (L,W,H) and centered (xc,yc,zc), which
> is much smaller than the model domain, and we set Vx = Vpush, Vy=0 within
> the block (Vz is let free for easier convergence).
> > As I said before, since the code does not have a monolithic matrix, but
> 4 submatrices (VV VP; PV PP), and the rhs has 2 sub vectors rhs=(f; g), my
> approach is to modify only (VV, VP, f) for the Dirichlet BC.
> >
> > The way I tested the implementation:
> > 1) Output (VV, VP, f, Dirichlet dofs) - unmodified (no Dirichlet BC)
> > 2) Output (VV, VP, f, Dirichlet dofs) - a) modified with MatZeroRowsIS(),
> > - b) modified with MatZeroRowsColumnsIS() -> S_PETSc
> > Again, the only difference between a) and b) is:
> > //
> > ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
> > // ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs);
> CHKERRQ(ierr);
> >
> > ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
> > ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
> > 3) Read them in Matlab and perform the exact same operations on the
> unmodified matrices and f vector. -> S_Matlab
> > 4) Compare S_PETSc with S_Matlab. If the implementation is correct, they
> should be equal (VV, VP, f).
> > 5) Check for 1 cpu and 4 cpus.
> >
> > Now to answer your questions:
> >
> > a,b,d) Yes, matrix modification is done correctly (check the spy
> diagrams below) in all cases: MatZeroRowsIS() and MatZeroRowsColumnsIS()
> on 1 and 4 cpus.
> >
> > I should have said that in the piece of code above:
> > v_vv = 1.0;
> > v_vp = 0.0;
> > The vector x_push is a duplicate of rhs, with zero elements except the
> values for the Dirichlet dofs.
> >
> > c) The rhs is a different matter. With MatZeroRows() there is no
> problem. The rhs is equivalent with the one in Matlab, sequential and
> parallel.
> > However, with MatZeroRowsColumns(), the residual contains nonzero
> elements, and in parallel the nonzero pattern is even bigger (1 cpu - 63, 4
> cpu - 554). But if you look carefully, the values of the nonzero residuals
> are very small < +/- 1e-10.
> > So, I did a tolerance filter:
> >
> > tol = 1e-10;
> > res = f_petsc - f_mod_matlab;
> > for i=1:length(res)
> > if abs(res(i))>0 & abs(res(i))<tol
> > res(i)=0;
> > end
> > end
> >
> > and then the f_petsc and f_mod_matlab are equivalent on 1 and 4 cpus
> (figure 5). So it seems that MatZeroRowsColumnsIS() might give some nonzero
> residuals.
> >
> > Calculating the norm2 of the residuals defined above in each case gives:
> > MatZeroRowsIS() 1cpu: norm(res,2) = 0
> > MatZeroRowsIS() 4cpu: norm(res,2) = 0
> > MatZeroRowsColumnsIS() 1cpu: norm(res,2) = 1.6880e-10
> > MatZeroRowsColumnsIS() 4cpu: norm(res,2) = 7.3786e+06
> >
> > Since this is purely a problem of matrix and vector
> assembly/manipulation, I think the nonzero residuals of the rhs with
> MatZeroRowsColumnsIS() give the parallel artefacts that I showed last time.
> > If you need the raw data and the matlab scripts that I used for testing
> for your consideration, please let me know.
> >
> > Thanks,
> > Adina
> >
> > When performing the manual operations on the unmodified matrices and rhs
> vector in Matlab, I took into account:
> > - matlab indexing = petsc indexing +1;
> > - the vectors written to file for matlab (PETSC_VIEWER_BINARY_MATLAB)
> have the natural ordering, rather than the petsc ordering. On 1 cpu, they
> are equivalent, but on 4 cpus, the Dirichlet BC indices had to be converted
> to natural indices in order to perform the correct operations on the rhs.
> >
> > <spy_zerorows_1cpu.png>
> > <spy_zerorows_4cpu.png>
> > <spy_zerorowscol_1cpu.png>
> > <spy_zerorowscol_4cpu.png>
> > <residual_tol.png>
> >
> > On May 6, 2014, at 4:22 PM, Matthew Knepley wrote:
> >
> >> On Tue, May 6, 2014 at 7:23 AM, Püsök, Adina-Erika <
> puesoek at uni-mainz.de> wrote:
> >> Hello!
> >>
> >> I was trying to implement some internal Dirichlet boundary conditions
> into an aij matrix of the form: A=( VV VP; PV PP ). The idea was to
> create an internal block (let's say Dirichlet block) that moves with
> constant velocity within the domain (i.e. check all the dofs within the
> block and set the values accordingly to the desired motion).
> >>
> >> Ideally, this means to zero the rows and columns in VV, VP, PV
> corresponding to the dirichlet dofs and modify the corresponding rhs
> values. However, since we have submatrices and not a monolithic matrix A,
> we can choose to modify only VV and PV matrices.
> >> The global indices of the velocity points within the Dirichlet block
> are contained in the arrays rowid_array.
> >>
> >> What I want to point out is that the function MatZeroRowsColumnsIS()
> seems to create parallel artefacts, compared to MatZeroRowsIS() when run on
> more than 1 processor. Moreover, the results on 1 cpu are identical.
> >> See below the results of the test (the Dirichlet block is outlined in
> white) and the piece of the code involved where the 1) - 2) parts are the
> only difference.
> >>
> >> I am assuming that you are showing the result of solving the equations.
> It would be more useful, and presumably just as easy
> >> to say:
> >>
> >> a) Are the correct rows zeroed out?
> >>
> >> b) Is the diagonal element correct?
> >>
> >> c) Is the rhs value correct?
> >>
> >> d) Are the columns zeroed correctly?
> >>
> >> If we know where the problem is, its easier to fix. For example, if the
> rhs values are
> >> correct and the rows are zeroed, then something is wrong with the
> solution procedure.
> >> Since ZeroRows() works and ZeroRowsColumns() does not, this is a
> distinct possibility.
> >>
> >> Thanks,
> >>
> >> Matt
> >>
> >> Thanks,
> >> Adina Pusok
> >>
> >> // Create an IS required by MatZeroRows()
> >> ierr =
> ISCreateGeneral(PETSC_COMM_WORLD,numRowsx,rowidx_array,PETSC_COPY_VALUES,&isx);
> CHKERRQ(ierr);
> >> ierr =
> ISCreateGeneral(PETSC_COMM_WORLD,numRowsy,rowidy_array,PETSC_COPY_VALUES,&isy);
> CHKERRQ(ierr);
> >> ierr =
> ISCreateGeneral(PETSC_COMM_WORLD,numRowsz,rowidz_array,PETSC_COPY_VALUES,&isz);
> CHKERRQ(ierr);
> >>
> >> 1) /*
> >> ierr = MatZeroRowsColumnsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
> >> ierr = MatZeroRowsColumnsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
> >> ierr = MatZeroRowsColumnsIS(VV_MAT,isz,v_vv,x_push,rhs);
> CHKERRQ(ierr);*/
> >>
> >> 2) ierr = MatZeroRowsIS(VV_MAT,isx,v_vv,x_push,rhs); CHKERRQ(ierr);
> >> ierr = MatZeroRowsIS(VV_MAT,isy,v_vv,x_push,rhs); CHKERRQ(ierr);
> >> ierr = MatZeroRowsIS(VV_MAT,isz,v_vv,x_push,rhs); CHKERRQ(ierr);
> >>
> >> ierr = MatZeroRowsIS(VP_MAT,isx,v_vp,PETSC_NULL,PETSC_NULL);
> CHKERRQ(ierr);
> >> ierr = MatZeroRowsIS(VP_MAT,isy,v_vp,PETSC_NULL,PETSC_NULL);
> CHKERRQ(ierr);
> >> ierr = MatZeroRowsIS(VP_MAT,isz,v_vp,PETSC_NULL,PETSC_NULL);
> CHKERRQ(ierr);
> >>
> >> ierr = ISDestroy(&isx); CHKERRQ(ierr);
> >> ierr = ISDestroy(&isy); CHKERRQ(ierr);
> >> ierr = ISDestroy(&isz); CHKERRQ(ierr);
> >>
> >>
> >> Results (velocity) with MatZeroRowsColumnsIS().
> >> 1cpu<r01_1cpu_rows_columns.png> 4cpu<r01_rows_columns.png>
> >>
> >> Results (velocity) with MatZeroRowsIS():
> >> 1cpu<r01_1cpu_rows.png> 4cpu<r01_rows.png>
> >>
> >>
> >>
> >> --
> >> 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
> >
>
>
--
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/20140514/e2aba086/attachment.html>
More information about the petsc-users
mailing list