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

[Samuel Lanthaler]

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
    [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/2e1a824d/attachment.html>