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
>
