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
