[petsc-users] MatPtAP problem after subsequent call to MatDuplicate

[Hong]

Samuel:
It is fixed
https://bitbucket.org/petsc/petsc/commits/cff31925ac4fa731f75d96f7dbc9974207834be9

It will be merged to petsc-release once it passes all regression tests.
Thanks for your report!

Hong

Thank you for your swift reply, Hong!
>
> Yes, what you said is exactly what I'm doing. No, I don't get an error
> when doing a MatView of matC.
>
> Samuel
>
>
> On 12/05/2017 04:38 PM, Hong wrote:
>
> Samuel:
> You try to do following:
> 1) Create A;
> 2) Create P;
> 3) C = PtAP:
>      CALL MatPtAP(matA,matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,matC
> ,ierr);
> 4) MatDestroy(matA,ierr);
> 5) MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr);
> 6) MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr);
>
> The error occurs at (6). Do you get any error for
> MatView(matC)?
>
> C has some special data structures as a matrix product which could be lost
> during
> MatDuplicate for matA. It might be a bug in our code. I'll check it.
>
> Hong
>
>
>> Hi there,
>>
>> I am getting error messages after using MatPtAP to create a new matrix C
>> = Pt*A*P and then trying to assign A=C. The following is a minimal example
>> reproducing the problem:
>>
>> #include "slepc/finclude/slepc.h"
>>       USE slepcsys
>>       USE slepceps
>>       IMPLICIT NONE
>>       LOGICAL :: cause_error
>>       ! --- pure PETSc
>>       PetscErrorCode :: ierr
>>       PetscScalar :: vals(3,3), val
>>       Mat :: matA, matP, matC
>>       PetscInt :: m,idxn(3),idxm(3),idone(1)
>>
>>       ! initialize SLEPc & Petsc etc.
>>       CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
>>
>>       ! Set up a new matrix
>>       m = 3
>>       CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);
>>       CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);
>>       CALL MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);
>> CHKERRQ(ierr);
>>       CALL MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr);
>> CHKERRQ(ierr);
>>
>>       ! [.... call to MatSetValues to set values of matA]
>>
>>       ! assemble matrix
>>       CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
>>       CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
>>
>>       ! duplicate matrix
>>       CALL MatDuplicate(matA,MAT_DO_NOT_COPY_VALUES,matP,ierr);
>> CHKERRQ(ierr);
>>
>>      ! [.... call to MatSetValues to set values of matP]
>>
>>       ! assemble matrix
>>       CALL MatAssemblyBegin(matP,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
>>       CALL MatAssemblyEnd(matP,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr);
>>
>>       ! compute C=Pt*A*P
>>       cause_error = .TRUE. ! set to .TRUE. to cause error, .FALSE. seems
>> to work fine
>>       IF(.NOT.cause_error) THEN
>>          CALL MatDuplicate(matA,MAT_COPY_VALUES,matC,ierr);
>> CHKERRQ(ierr);
>>       ELSE
>>          CALL MatPtAP(matA,matP,MAT_INITIAL_
>> MATRIX,PETSC_DEFAULT_REAL,matC,ierr); CHKERRQ(ierr);
>>       END IF
>>
>>       ! destroy matA and duplicate A=C
>>       CALL MatDestroy(matA,ierr); CHKERRQ(ierr);
>>       CALL MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr); CHKERRQ(ierr);
>>
>>       ! display resulting matrix A
>>       CALL MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr); CHKERRQ(ierr);
>>
>> Whether A and P are assigned any values doesn't seem to matter at all.
>> The error message I'm getting is:
>>
>> Mat Object: 1 MPI processes
>>   type: mpiaij
>> [0]PETSC ERROR: ------------------------------
>> ------------------------------------------
>> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
>> probably memory access out of range
>> [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
>> [0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/d
>> ocumentation/faq.html#valgrind
>> [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS
>> X to find memory corruption errors
>> [0]PETSC ERROR: likely location of problem given in stack below
>> [0]PETSC ERROR: ---------------------  Stack Frames
>> ------------------------------------
>> [0]PETSC ERROR: Note: The EXACT line numbers in the stack are not
>> available,
>> [0]PETSC ERROR:       INSTEAD the line number of the start of the function
>> [0]PETSC ERROR:       is given.
>> [0]PETSC ERROR: [0] MatView_MPIAIJ_PtAP line 23
>> /home/lanthale/Progs/petsc-3.8.2/src/mat/impls/aij/mpi/mpiptap.c
>> [0]PETSC ERROR: [0] MatView line 949 /home/lanthale/Progs/petsc-3.8
>> .2/src/mat/interface/matrix.c
>> [0]PETSC ERROR: --------------------- Error Message
>> --------------------------------------------------------------
>> [0]PETSC ERROR: Signal received
>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>> for trouble shooting.
>> [0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017
>>
>> Could someone maybe tell me where I'm doing things wrong? The documentation
>> for MatPtAP
>> <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAP.html>
>> says that C will be created and that this routine is "currently only
>> implemented for pairs of AIJ matrices and classes which inherit from AIJ".
>> Does this maybe exclude MPIAIJ matrices?
>>
>> Thanks,
>> Samuel
>>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171205/4c8900aa/attachment-0001.html>