<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Samuel:</div><div class="gmail_quote">You try to do following:</div><div class="gmail_quote">1) Create A;</div><div class="gmail_quote">2) Create P;</div><div class="gmail_quote">3) C = PtAP:</div><div class="gmail_quote"> <span style="font-size:10.6667px">CALL MatPtAP(matA,matP,MAT_INITIAL_</span><wbr style="font-size:10.6667px"><span style="font-size:10.6667px">MATRIX,PETSC_DEFAULT_REAL,</span><wbr style="font-size:10.6667px"><span style="font-size:10.6667px">matC,ierr); </span></div><div class="gmail_quote"><span style="font-size:10.6667px">4) </span><span style="font-size:10.6667px">MatDestroy(matA,ierr);</span></div><div class="gmail_quote"><span style="font-size:10.6667px">5) </span><span style="font-size:10.6667px">MatDuplicate(matC,MAT_COPY_</span><wbr style="font-size:10.6667px"><span style="font-size:10.6667px">VALUES,matA,ierr); </span></div><div class="gmail_quote"><span style="font-size:10.6667px">6) </span><span style="font-size:10.6667px">MatView(matA,PETSC_VIEWER_</span><wbr style="font-size:10.6667px"><span style="font-size:10.6667px">STDOUT_WORLD,ierr); </span></div><div class="gmail_quote"><span style="font-size:10.6667px"><br></span></div><div class="gmail_quote"><span style="font-size:10.6667px">The error occurs at (6). Do you get any error for</span></div><div class="gmail_quote"><span style="font-size:10.6667px">MatView(matC)?</span></div><div class="gmail_quote"><br></div><div class="gmail_quote">C has some special data structures as a matrix product which could be lost during</div><div class="gmail_quote">MatDuplicate for matA. It might be a bug in our code. I'll check it.</div><div class="gmail_quote"><br></div><div class="gmail_quote">Hong</div><div class="gmail_quote"><span style="font-size:10.6667px"> <br></span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div bgcolor="#FFFFFF">
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_<wbr>CHARACTER,ierr)<br>
<br>
! Set up a new matrix<br>
m = 3<br>
CALL MatCreate(PETSC_COMM_WORLD,<wbr>matA,ierr); CHKERRQ(ierr);<br>
CALL MatSetType(matA,MATMPIAIJ,<wbr>ierr); CHKERRQ(ierr);<br>
CALL MatSetSizes(matA,PETSC_DECIDE,<wbr>PETSC_DECIDE,m,m,ierr);
CHKERRQ(ierr);<br>
CALL
MatMPIAIJSetPreallocation(<wbr>matA,3,PETSC_NULL_INTEGER,3,<wbr>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_<wbr>FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
CALL MatAssemblyEnd(matA,MAT_FINAL_<wbr>ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
<br>
! duplicate matrix<br>
CALL MatDuplicate(matA,MAT_DO_NOT_<wbr>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_<wbr>FINAL_ASSEMBLY,ierr);
CHKERRQ(ierr);<br>
CALL MatAssemblyEnd(matP,MAT_FINAL_<wbr>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_<wbr>VALUES,matC,ierr);
CHKERRQ(ierr);<br>
ELSE<br>
CALL
MatPtAP(matA,matP,MAT_INITIAL_<wbr>MATRIX,PETSC_DEFAULT_REAL,<wbr>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_<wbr>VALUES,matA,ierr);
CHKERRQ(ierr);<br>
<br>
! display resulting matrix A<br>
CALL MatView(matA,PETSC_VIEWER_<wbr>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:
------------------------------<wbr>------------------------------<wbr>------------<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="gmail-m_349307368681188326moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>documentation/faq.html#<wbr>valgrind</a><br>
[0]PETSC ERROR: or try <a class="gmail-m_349307368681188326moz-txt-link-freetext" href="http://valgrind.org" target="_blank">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
------------------------------<wbr>------<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.<wbr>8.2/src/mat/impls/aij/mpi/<wbr>mpiptap.c<br>
[0]PETSC ERROR: [0] MatView line 949
/home/lanthale/Progs/petsc-3.<wbr>8.2/src/mat/interface/matrix.c<br>
[0]PETSC ERROR: --------------------- Error Message
------------------------------<wbr>------------------------------<wbr>--<br>
[0]PETSC ERROR: Signal received<br>
[0]PETSC ERROR: See
<a class="gmail-m_349307368681188326moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>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" target="_blank">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>
</div>
</blockquote></div><br></div></div>