Pierre : I'm not aware the "dummy function" MatMatMultNumeric_MPIDense(). Good to know you find a solution.
Hong > Barry, > There is no such error when you do: > --- a/src/mat/examples/tests/ex109.c > +++ b/src/mat/examples/tests/ex109.c > @@ -57,7 +57,8 @@ int main(int argc,char **argv) > ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); > > /* Test C = A*B (aij*dense) */ > - ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); > + // ierr = MatMatMult(A,B,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); > + MatCreateDense(PETSC_COMM_WORLD,an,PETSC_DECIDE,PETSC_ > DECIDE,M,NULL,&C); > > What you are showing is maybe an error in the “MatDuplicate” that > duplicates the ops->matmultnumeric, but not the workB container? > > Hong, > Once again, I’ve been using MatMatMult without any calls to > MatMatMultSymbolic_MPIBAIJ_MPIDense for many months. > What you are saying is only partially true, the workspace array workB may > be allocated in MatMatMultSymbolic_MPIBAIJ_MPIDense, but also in the > routine MatMatMultNumeric_MPIDense which exact role, as written in the > source code, is to bypass MatMatMultSymbolic_MPIBAIJ_MPIDense: > /* > This is a "dummy function" that handles the case where matrix C was > created as a dense matrix > directly by the user and passed to MatMatMult() with the > MAT_REUSE_MATRIX option > > It is the same as MatMatMultSymbolic_MPIAIJ_MPIDense() except does not > create C > */ > > Anyway, I realised that my question was dumb in the first place. In > MatMatMultNumeric_MPIDense, I just check now when A is MATMPIBAIJ or > MATMPIAIJ and I either cast A to Mat_MPIBAIJ* or Mat_MPIAIJ* and set > C->ops->matmultnumeric accordingly. > Now, everything works just fine, with *no* prior call to MatMatMult with > MAT_INITIAL_MATRIX. > > Thanks for your patience, > Pierre > > > On 8 Jun 2017, at 5:14 AM, Zhang, Hong <[email protected]> wrote: > > > > Pierre: > > You must call MatMatMultSymbolic_MPIBAIJ_MPIDense() to create a dense > matrix C (=A*B) and an internal submatrix 'workB'. > > 'workB' will be used (and reused) to store off processor rows of B > needed for local product of C. > > See MatMatMultSymbolic_MPIAIJ_MPIDense(). > > > > Without this local storage, numeric product fails as shown from your > test. > > Again, for using parallel MatMatMult(), you must use > MatMatMultSymbolic() to create C, then reuse it. > > > > Hong > > ________________________________________ > > From: Barry Smith [[email protected]] > > Sent: Wednesday, June 07, 2017 9:31 PM > > To: Pierre Jolivet > > Cc: Zhang, Hong; For users of the development version of PETSc > > Subject: Re: [petsc-dev] MatMatMult > > > >> On Jun 7, 2017, at 1:29 AM, Pierre Jolivet <[email protected]> > wrote: > >> > >>> > >>> On 6 Jun 2017, at 11:21 PM, Barry Smith <[email protected]> wrote: > >>> > >>>> > >>>> On Jun 6, 2017, at 1:52 PM, Pierre Jolivet < > [email protected]> wrote: > >>>> > >>>> On Tue, 6 Jun 2017 13:44:06 -0500, Hong wrote: > >>>>> Pierre : > >>>>> > >>>>>> Yes, of course I defined (*C)->ops->matmultnumeric = > >>>>>> MatMatMultNumeric_MPIBAIJ_MPIDense in > >>>>>> MatMatMultSymbolic_MPIBAIJ_MPIDense. > >>>>>> However, the routine MatMatMultSymbolic_MPIBAIJ_MPIDense is never > >>>>>> reached when calling MatMatMult with scall == MAT_REUSE_MATRIX > >>>>>> > >>>>> (http://www.mcs.anl.gov/petsc/petsc-current/src/mat/ > interface/matrix.c.html#line9487 > >>>>>> [1], just to be sure I added a dummy printf in > >>>>>> MatMatMultSymbolic_MPIBAIJ_MPIDense and nothing is displayed with > >>>>>> MAT_REUSE_MATRIX, MAT_INITIAL_MATRIX works as intended) > >>>>> > >>>>> MatMatMultSymbolic_xxx() is called only for MAT_INITIAL_MATRIX, > >>>>> during which, it defines > >>>>> (*C)->ops->matmultnumeric = MatMatMultNumeric_MPIBAIJ_MPIDense; > >>>> > >>>> Sorry, I was not thorough enough. What I meant to say is that I never > call MatMatMult with MAT_INITIAL_MATRIX, so MatMatMultSymbolic_xxx is never > called. > >>> > >>> It is not designed to work that way. The reuse is specifically to > reuse a matrix obtained with a call to initial_matrix, not to reuse a > random matrix the user provides. The reason is that the symbolic stage > may/does build additional data structures that are used for the multiples > that the user cannot see or build directly. > >> > >> I’m doing: > >> MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, in, &B); > >> MatCreateDense(PETSC_COMM_WORLD, m, n, M, N, out, &C); > >> MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C); > >> > >> This works perfectly fine with A of type MPIAIJ, but once again, > everything is hardwired for this only type. > > > > Hmm, I modified src/mat/examples/tests/ex109.c to create the dense D > matrix with a MatDuplicate() instead of INITIAL MATRIX and when I run in > parallel I get > > > > $ petscmpiexec -n 2 ./ex109 > > [0]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > > [0]PETSC ERROR: Petsc has generated inconsistent data > > [0]PETSC ERROR: Container does not exist > > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > > [0]PETSC ERROR: Petsc Release Version 3.7.6, unknown > > [0]PETSC ERROR: ./ex109 on a arch-basic named Barrys-MacBook-Pro.local > by barrysmith Wed Jun 7 21:24:59 2017 > > [0]PETSC ERROR: Configure options PETSC_ARCH=arch-basic > > [0]PETSC ERROR: #1 MatMPIDenseScatter() line 487 in > /Users/barrysmith/Src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c > > [0]PETSC ERROR: #2 MatMatMultNumeric_MPIAIJ_MPIDense() line 557 in > /Users/barrysmith/Src/petsc/src/mat/impls/aij/mpi/mpimatmatmult.c > > [0]PETSC ERROR: #3 MatMatMult() line 9492 in /Users/barrysmith/Src/petsc/ > src/mat/interface/matrix.c > > [0]PETSC ERROR: #4 main() line 79 in /Users/barrysmith/Src/petsc/ > src/mat/examples/tests/ex109.c > > [0]PETSC ERROR: PETSc Option Table entries: > > [0]PETSC ERROR: -malloc_test > > > > it does run sequentially. > > > >> According to your answer, Barry, I think the manual should be adjusted > accordingly because to me what I’m trying to achieve should be possible as > is (the fact that it is possible with MPIAIJ is also puzzling): > >> "In the special case where matrix B (and hence C) are dense you can > create the correctly sized matrix C yourself and then call this routine > with MAT_REUSE_MATRIX, rather than first having MatMatMult() create it for > you. […]” (http://www.mcs.anl.gov/petsc/petsc-current/docs/ > manualpages/Mat/MatMatMult.html) > > > > Yes, you are right this documentation is wrong. This is only true for > sequential dense matrices. I've fixed the manual page in master. > > > > Barry > > > >> > >>>> Is there another way? > >>>> I cannot use MAT_INITIAL_MATRIX because I want to manage the > underlying memory, and I'm guessing MAT_INITIAL_MATRIX will delete whatever > is in C and reallocate memory on its own. > >>> > >>> The memory in C is just a dense array so presumably you want to use > some "application" memory for this space? > >>> > >>> The simplest way for you to achieve this is to call the > MatMatMultSymbolic() then dig out the underly array of the matrix, call > PetscFree() on it and put in your own array. Yes, this is slightly hacky, > but it will work fine within the current PETSc model. > >> > >> Yeah, alright, I already need to tinker with the > >> MatDensePlaceArray/MatDenseResetArray > so I can figure something out as well for this. > >> Thanks, > >> Pierre > >> > >>> Barry > >>> > >>>> > >>>> Thanks in advance for your help, > >>>> Pierre > >>>> > >>>>> Then MatMatMult(A,C,MAT_REUSE_MATRIX,..) calls > >>>>> (*(*C)->ops->matmultnumeric)(A,B,*C); (line 9432 in matrix.c) > >>>>> > >>>>> which should go to MatMatMultNumeric_MPIBAIJ_MPIDense. > >>>>> > >>>>> You may follow a debugging process > >>>>> using petsc/src/mat/examples/tests/ex109.c > >>>>> > >>>>> Are you working on a branch of petsc? If so, I may take a look and > >>>>> see what is the problem. > >>>>> > >>>>> Hong > >>>>> > >>>>>>>> 2) I'm having trouble when scall == MAT_REUSE_MATRIX. Here, > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>> > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > dense/mpi/mpidense.c.html#line1208 > >>>>>>> [2] > >>>>>>> > >>>>>>>> [2] it looks that the numeric part of the MatMatMult (which is > >>>>>>>> called when scall == MAT_REUSE_MATRIX) is hardwired to this > >>>>>>>> routine > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>> > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > aij/mpi/mpimatmatmult.c.html#line376 > >>>>>>> [3] > >>>>>>> [3]. > >>>>>>> br> > >>>>>>> > >>>>>> > >>>>> s.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/ > mpimatmatmult.c.html#line556" > >>>>>>> rel="noreferrer" > >>>>>>> > >>>>>> > >>>>> target="_blank">http://www.mcs.anl.gov/petsc/petsc- > current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line556 > >>>>>>> [2] > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>> > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > dense/mpi/mpidense.c.html#line1208 > >>>>>>> [4] > >>>>>>> [3] > >>>>>>> > >>>>>>> > >>>>>>> > >>>>>> > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > aij/mpi/mpimatmatmult.c.html#line376 > >>>>>>> [5] > >>>>> > >>>>> > >>>>> > >>>>> Links: > >>>>> ------ > >>>>> [1] > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/ > interface/matrix.c.html#line9487 > >>>>> [2] > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > dense/mpi/mpidense.c.html#line1208 > >>>>> [3] > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > aij/mpi/mpimatmatmult.c.html#line376 > >>>>> [4] > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > dense/mpi/mpidense.c.html#line1208 > >>>>> [5] > >>>>> http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/ > aij/mpi/mpimatmatmult.c.html#line376 > > > >
