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

[Samuel Lanthaler]

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/documentation/faq.html#valgrind
>         <http://www.mcs.anl.gov/petsc/documentation/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
>         <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/5961a53d/attachment.html>