<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=utf-8">
</head>
<body bgcolor="#FFFFFF" text="#000000">
Hi there,<br>
<br>
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:<br>
<blockquote><small>#include "slepc/finclude/slepc.h"<br>
USE slepcsys<br>
USE slepceps<br>
IMPLICIT NONE<br>
LOGICAL :: cause_error<br>
! --- pure PETSc<br>
PetscErrorCode :: ierr<br>
PetscScalar :: vals(3,3), val<br>
Mat :: matA, matP, matC<br>
PetscInt :: m,idxn(3),idxm(3),idone(1)<br>
<br>
! initialize SLEPc & Petsc etc.<br>
CALL SlepcInitialize(PETSC_NULL_CHARACTER,ierr)<br>
<br>
! Set up a new matrix<br>
m = 3<br>
CALL MatCreate(PETSC_COMM_WORLD,matA,ierr); CHKERRQ(ierr);<br>
CALL MatSetType(matA,MATMPIAIJ,ierr); CHKERRQ(ierr);<br>
CALL MatSetSizes(matA,PETSC_DECIDE,PETSC_DECIDE,m,m,ierr);
CHKERRQ(ierr);<br>
CALL
MatMPIAIJSetPreallocation(matA,3,PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,ierr);
CHKERRQ(ierr);<br>
<br>
! [.... call to MatSetValues to set values of matA]<br>
<br>
! assemble matrix<br>
CALL MatAssemblyBegin(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
CALL MatAssemblyEnd(matA,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
<br>
! duplicate matrix<br>
CALL MatDuplicate(matA,MAT_DO_NOT_COPY_VALUES,matP,ierr);
CHKERRQ(ierr);<br>
<br>
! [.... call to MatSetValues to set values of matP]<br>
<br>
! assemble matrix<br>
CALL MatAssemblyBegin(matP,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
CALL MatAssemblyEnd(matP,MAT_FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
<br>
! compute C=Pt*A*P<br>
cause_error = .TRUE. ! set to .TRUE. to cause error,
.FALSE. seems to work fine<br>
IF(.NOT.cause_error) THEN<br>
CALL MatDuplicate(matA,MAT_COPY_VALUES,matC,ierr);
CHKERRQ(ierr);<br>
ELSE<br>
CALL
MatPtAP(matA,matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,matC,ierr);
CHKERRQ(ierr);<br>
END IF<br>
<br>
! destroy matA and duplicate A=C<br>
CALL MatDestroy(matA,ierr); CHKERRQ(ierr);<br>
CALL MatDuplicate(matC,MAT_COPY_VALUES,matA,ierr);
CHKERRQ(ierr);<br>
<br>
! display resulting matrix A<br>
CALL MatView(matA,PETSC_VIEWER_STDOUT_WORLD,ierr);
CHKERRQ(ierr);</small><br>
<br>
</blockquote>
Whether A and P are assigned any values doesn't seem to matter at
all. The error message I'm getting is:<br>
<blockquote><small>Mat Object: 1 MPI processes<br>
type: mpiaij<br>
[0]PETSC ERROR:
------------------------------------------------------------------------<br>
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation
Violation, probably memory access out of range<br>
[0]PETSC ERROR: Try option -start_in_debugger or
-on_error_attach_debugger<br>
[0]PETSC ERROR: or see
<a class="moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind</a><br>
[0]PETSC ERROR: or try <a class="moz-txt-link-freetext" href="http://valgrind.org">http://valgrind.org</a> on GNU/linux and
Apple Mac OS X to find memory corruption errors<br>
[0]PETSC ERROR: likely location of problem given in stack below<br>
[0]PETSC ERROR: --------------------- Stack Frames
------------------------------------<br>
[0]PETSC ERROR: Note: The EXACT line numbers in the stack are
not available,<br>
[0]PETSC ERROR: INSTEAD the line number of the start of
the function<br>
[0]PETSC ERROR: is given.<br>
[0]PETSC ERROR: [0] MatView_MPIAIJ_PtAP line 23
/home/lanthale/Progs/petsc-3.8.2/src/mat/impls/aij/mpi/mpiptap.c<br>
[0]PETSC ERROR: [0] MatView line 949
/home/lanthale/Progs/petsc-3.8.2/src/mat/interface/matrix.c<br>
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------<br>
[0]PETSC ERROR: Signal received<br>
[0]PETSC ERROR: See
<a class="moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble
shooting.<br>
[0]PETSC ERROR: Petsc Release Version 3.8.2, Nov, 09, 2017 </small><br>
<br>
</blockquote>
Could someone maybe tell me where I'm doing things wrong? The <a
href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAP.html">documentation
for MatPtAP</a> 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?<br>
<br>
Thanks,<br>
Samuel<br>
</body>
</html>