<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Pierre :</div><div class="gmail_quote">I'm not aware the "dummy function" MatMatMultNumeric_MPIDense(). </div><div class="gmail_quote">Good to know you find a solution.</div><div class="gmail_quote"><br></div><div class="gmail_quote">Hong<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Barry,<br>
There is no such error when you do:<br>
--- a/src/mat/examples/tests/<wbr>ex109.c<br>
+++ b/src/mat/examples/tests/<wbr>ex109.c<br>
@@ -57,7 +57,8 @@ int main(int argc,char **argv)<br>
   ierr = MatAssemblyEnd(B,MAT_FINAL_<wbr>ASSEMBLY);CHKERRQ(ierr);<br>
<br>
   /* Test C = A*B (aij*dense) */<br>
-  ierr = MatMatMult(A,B,MAT_INITIAL_<wbr>MATRIX,fill,&C);CHKERRQ(ierr);<br>
+  // ierr = MatMatMult(A,B,MAT_INITIAL_<wbr>MATRIX,fill,&C);CHKERRQ(ierr);<br>
+  MatCreateDense(PETSC_COMM_<wbr>WORLD,an,PETSC_DECIDE,PETSC_<wbr>DECIDE,M,NULL,&C);<br>
<br>
What you are showing is maybe an error in the “MatDuplicate” that duplicates the ops->matmultnumeric, but not the workB container?<br>
<br>
Hong,<br>
Once again, I’ve been using MatMatMult without any calls to MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense for many months.<br>
What you are saying is only partially true, the workspace array workB may be allocated in MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense, but also in the routine MatMatMultNumeric_MPIDense which exact role, as written in the source code, is to bypass MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense:<br>
/*<br>
    This is a "dummy function" that handles the case where matrix C was created as a dense matrix<br>
  directly by the user and passed to MatMatMult() with the MAT_REUSE_MATRIX option<br>
<br>
  It is the same as MatMatMultSymbolic_MPIAIJ_<wbr>MPIDense() except does not create C<br>
*/<br>
<br>
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.<br>
Now, everything works just fine, with *no* prior call to MatMatMult with MAT_INITIAL_MATRIX.<br>
<br>
Thanks for your patience,<br>
Pierre<br>
<div class="gmail-HOEnZb"><div class="gmail-h5"><br>
> On 8 Jun 2017, at 5:14 AM, Zhang, Hong <<a href="mailto:hzhang@mcs.anl.gov">hzhang@mcs.anl.gov</a>> wrote:<br>
><br>
> Pierre:<br>
> You must call MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense() to create a dense matrix C (=A*B) and an internal submatrix 'workB'.<br>
> 'workB' will be used (and reused) to store off processor rows of B needed for local product of C.<br>
> See MatMatMultSymbolic_MPIAIJ_<wbr>MPIDense().<br>
><br>
> Without this local storage, numeric product fails as shown from your test.<br>
> Again, for using parallel MatMatMult(), you must use MatMatMultSymbolic() to create C, then reuse it.<br>
><br>
> Hong<br>
> ______________________________<wbr>__________<br>
> From: Barry Smith [<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>]<br>
> Sent: Wednesday, June 07, 2017 9:31 PM<br>
> To: Pierre Jolivet<br>
> Cc: Zhang, Hong; For users of the development version of PETSc<br>
> Subject: Re: [petsc-dev] MatMatMult<br>
><br>
>> On Jun 7, 2017, at 1:29 AM, Pierre Jolivet <<a href="mailto:pierre.jolivet@enseeiht.fr">pierre.jolivet@enseeiht.fr</a>> wrote:<br>
>><br>
>>><br>
>>> On 6 Jun 2017, at 11:21 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>>><br>
>>>><br>
>>>> On Jun 6, 2017, at 1:52 PM, Pierre Jolivet <<a href="mailto:Pierre.Jolivet@enseeiht.fr">Pierre.Jolivet@enseeiht.fr</a>> wrote:<br>
>>>><br>
>>>> On Tue, 6 Jun 2017 13:44:06 -0500, Hong wrote:<br>
>>>>> Pierre :<br>
>>>>><br>
>>>>>> Yes, of course I defined (*C)->ops->matmultnumeric =<br>
>>>>>> MatMatMultNumeric_MPIBAIJ_<wbr>MPIDense in<br>
>>>>>> MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense.<br>
>>>>>> However, the routine MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense is never<br>
>>>>>> reached when calling MatMatMult with scall == MAT_REUSE_MATRIX<br>
>>>>>><br>
>>>>> (<a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface/matrix.c.html#line9487" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/<wbr>interface/matrix.c.html#<wbr>line9487</a><br>
>>>>>> [1], just to be sure I added a dummy printf in<br>
>>>>>> MatMatMultSymbolic_MPIBAIJ_<wbr>MPIDense and nothing is displayed with<br>
>>>>>> MAT_REUSE_MATRIX, MAT_INITIAL_MATRIX works as intended)<br>
>>>>><br>
>>>>> MatMatMultSymbolic_xxx() is called only for MAT_INITIAL_MATRIX,<br>
>>>>> during which, it defines<br>
>>>>> (*C)->ops->matmultnumeric = MatMatMultNumeric_MPIBAIJ_<wbr>MPIDense;<br>
>>>><br>
>>>> 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.<br>
>>><br>
>>>  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.<br>
>><br>
>> I’m doing:<br>
>>  MatCreateDense(PETSC_COMM_<wbr>WORLD, m, n, M, N, in, &B);<br>
>>  MatCreateDense(PETSC_COMM_<wbr>WORLD, m, n, M, N, out, &C);<br>
>>  MatMatMult(A, B, MAT_REUSE_MATRIX, PETSC_DEFAULT, &C);<br>
>><br>
>> This works perfectly fine with A of type MPIAIJ, but once again, everything is hardwired for this only type.<br>
><br>
>   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<br>
><br>
> $ petscmpiexec -n 2 ./ex109<br>
> [0]PETSC ERROR: --------------------- Error Message ------------------------------<wbr>------------------------------<wbr>--<br>
> [0]PETSC ERROR: Petsc has generated inconsistent data<br>
> [0]PETSC ERROR: Container does not exist<br>
> [0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>documentation/faq.html</a> for trouble shooting.<br>
> [0]PETSC ERROR: Petsc Release Version 3.7.6, unknown<br>
> [0]PETSC ERROR: ./ex109 on a arch-basic named Barrys-MacBook-Pro.local by barrysmith Wed Jun  7 21:24:59 2017<br>
> [0]PETSC ERROR: Configure options PETSC_ARCH=arch-basic<br>
> [0]PETSC ERROR: #1 MatMPIDenseScatter() line 487 in /Users/barrysmith/Src/petsc/<wbr>src/mat/impls/aij/mpi/<wbr>mpimatmatmult.c<br>
> [0]PETSC ERROR: #2 MatMatMultNumeric_MPIAIJ_<wbr>MPIDense() line 557 in /Users/barrysmith/Src/petsc/<wbr>src/mat/impls/aij/mpi/<wbr>mpimatmatmult.c<br>
> [0]PETSC ERROR: #3 MatMatMult() line 9492 in /Users/barrysmith/Src/petsc/<wbr>src/mat/interface/matrix.c<br>
> [0]PETSC ERROR: #4 main() line 79 in /Users/barrysmith/Src/petsc/<wbr>src/mat/examples/tests/ex109.c<br>
> [0]PETSC ERROR: PETSc Option Table entries:<br>
> [0]PETSC ERROR: -malloc_test<br>
><br>
> it does run sequentially.<br>
><br>
>> 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):<br>
>> "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. […]” (<a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMatMult.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/docs/<wbr>manualpages/Mat/MatMatMult.<wbr>html</a>)<br>
><br>
>   Yes, you are right this documentation is wrong. This is only true for sequential dense matrices. I've fixed the manual page in master.<br>
><br>
>  Barry<br>
><br>
>><br>
>>>> Is there another way?<br>
>>>> 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.<br>
>>><br>
>>> The memory in C is just a dense array so presumably you want to use some "application" memory for this space?<br>
>>><br>
>>>  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.<br>
>><br>
>> Yeah, alright, I already need to tinker with the MatDensePlaceArray/<wbr>MatDenseResetArray so I can figure something out as well for this.<br>
>> Thanks,<br>
>> Pierre<br>
>><br>
>>> Barry<br>
>>><br>
>>>><br>
>>>> Thanks in advance for your help,<br>
>>>> Pierre<br>
>>>><br>
>>>>> Then MatMatMult(A,C,MAT_REUSE_<wbr>MATRIX,..) calls<br>
>>>>> (*(*C)->ops->matmultnumeric)(<wbr>A,B,*C); (line 9432 in matrix.c)<br>
>>>>><br>
>>>>> which should go to MatMatMultNumeric_MPIBAIJ_<wbr>MPIDense.<br>
>>>>><br>
>>>>> You may follow a debugging process<br>
>>>>> using petsc/src/mat/examples/tests/<wbr>ex109.c<br>
>>>>><br>
>>>>> Are you working on a branch of petsc? If so, I may take a look and<br>
>>>>> see what is the problem.<br>
>>>>><br>
>>>>> Hong<br>
>>>>><br>
>>>>>>>> 2) I'm having trouble when scall == MAT_REUSE_MATRIX. Here,<br>
>>>>>>><br>
>>>>>>><br>
>>>>>>><br>
>>>>>><br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/dense/mpi/mpidense.c.html#line1208" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>dense/mpi/mpidense.c.html#<wbr>line1208</a><br>
>>>>>>> [2]<br>
>>>>>>><br>
>>>>>>>> [2] it looks that the numeric part of the MatMatMult (which is<br>
>>>>>>>> called when scall == MAT_REUSE_MATRIX) is hardwired to this<br>
>>>>>>>> routine<br>
>>>>>>><br>
>>>>>>><br>
>>>>>>><br>
>>>>>><br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line376" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>aij/mpi/mpimatmatmult.c.html#<wbr>line376</a><br>
>>>>>>> [3]<br>
>>>>>>> [3].<br>
>>>>>>> br><br>
>>>>>>><br>
>>>>>><br>
>>>>> <a href="http://s.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line556" rel="noreferrer" target="_blank">s.anl.gov/petsc/petsc-current/<wbr>src/mat/impls/aij/mpi/<wbr>mpimatmatmult.c.html#line556</a>"<br>
>>>>>>> rel="noreferrer"<br>
>>>>>>><br>
>>>>>><br>
>>>>> target="_blank"><a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line556" rel="noreferrer" target="_blank">http://www.<wbr>mcs.anl.gov/petsc/petsc-<wbr>current/src/mat/impls/aij/mpi/<wbr>mpimatmatmult.c.html#line556</a><br>
>>>>>>> [2]<br>
>>>>>>><br>
>>>>>>><br>
>>>>>>><br>
>>>>>><br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/dense/mpi/mpidense.c.html#line1208" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>dense/mpi/mpidense.c.html#<wbr>line1208</a><br>
>>>>>>> [4]<br>
>>>>>>> [3]<br>
>>>>>>><br>
>>>>>>><br>
>>>>>>><br>
>>>>>><br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line376" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>aij/mpi/mpimatmatmult.c.html#<wbr>line376</a><br>
>>>>>>> [5]<br>
>>>>><br>
>>>>><br>
>>>>><br>
>>>>> Links:<br>
>>>>> ------<br>
>>>>> [1]<br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/interface/matrix.c.html#line9487" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/<wbr>interface/matrix.c.html#<wbr>line9487</a><br>
>>>>> [2]<br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/dense/mpi/mpidense.c.html#line1208" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>dense/mpi/mpidense.c.html#<wbr>line1208</a><br>
>>>>> [3]<br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line376" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>aij/mpi/mpimatmatmult.c.html#<wbr>line376</a><br>
>>>>> [4]<br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/dense/mpi/mpidense.c.html#line1208" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>dense/mpi/mpidense.c.html#<wbr>line1208</a><br>
>>>>> [5]<br>
>>>>> <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/mat/impls/aij/mpi/mpimatmatmult.c.html#line376" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-current/src/mat/impls/<wbr>aij/mpi/mpimatmatmult.c.html#<wbr>line376</a><br>
><br>
<br>
</div></div></blockquote></div><br></div></div>