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 >> > > >
