<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Samuel:</div><div class="gmail_quote">It is fixed </div><div class="gmail_quote"><a href="https://bitbucket.org/petsc/petsc/commits/cff31925ac4fa731f75d96f7dbc9974207834be9">https://bitbucket.org/petsc/petsc/commits/cff31925ac4fa731f75d96f7dbc9974207834be9</a></div><div class="gmail_quote"><br></div><div class="gmail_quote">It will be merged to petsc-release once it passes all regression tests.</div><div class="gmail_quote">Thanks for your report!</div><div class="gmail_quote"><br></div><div class="gmail_quote">Hong</div><div class="gmail_quote"><br><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">
Thank you for your swift reply, Hong! <br>
<br>
Yes, what you said is exactly what I'm doing. No, I don't get an
error when doing a MatView of matC.<span class="gmail-HOEnZb"><font color="#888888"><br>
<br>
Samuel</font></span><div><div class="gmail-h5"><br>
<br>
<div class="gmail-m_-1435807948543973791moz-cite-prefix">On 12/05/2017 04:38 PM, Hong wrote:<br>
</div>
<blockquote type="cite">
<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><span style="font-size:10.6667px"><wbr>MATRIX,PETSC_DEFAULT_REAL,</span><span style="font-size:10.6667px">matC<wbr>,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><span style="font-size:10.6667px"><wbr>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><span style="font-size:10.6667px">S<wbr>TDOUT_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_CHA<wbr>RACTER,ierr)<br>
<br>
! Set up a new matrix<br>
m = 3<br>
CALL MatCreate(PETSC_COMM_WORLD,mat<wbr>A,ierr);
CHKERRQ(ierr);<br>
CALL MatSetType(matA,MATMPIAIJ,ierr<wbr>);
CHKERRQ(ierr);<br>
CALL MatSetSizes(matA,PETSC_DECIDE,<wbr>PETSC_DECIDE,m,m,ierr);
CHKERRQ(ierr);<br>
CALL MatMPIAIJSetPreallocation(matA<wbr>,3,PETSC_NULL_INTEGER,3,PETSC_<wbr>NULL_INTEGER,ierr);
CHKERRQ(ierr);<br>
<br>
! [.... call to MatSetValues to set values of
matA]<br>
<br>
! assemble matrix<br>
CALL MatAssemblyBegin(matA,MAT_FINA<wbr>L_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_C<wbr>OPY_VALUES,matP,ierr);
CHKERRQ(ierr);<br>
<br>
! [.... call to MatSetValues to set values of
matP]<br>
<br>
! assemble matrix<br>
CALL MatAssemblyBegin(matP,MAT_FINA<wbr>L_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_VAL<wbr>UES,matC,ierr);
CHKERRQ(ierr);<br>
ELSE<br>
CALL MatPtAP(matA,matP,MAT_INITIAL_<wbr>MATRIX,PETSC_DEFAULT_REAL,matC<wbr>,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_VAL<wbr>UES,matA,ierr);
CHKERRQ(ierr);<br>
<br>
! display resulting matrix A<br>
CALL MatView(matA,PETSC_VIEWER_STDO<wbr>UT_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_-1435807948543973791gmail-m_349307368681188326moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind" target="_blank">http://www.mcs.anl.gov/petsc/d<wbr>ocumentation/faq.html#valgrind</a><br>
[0]PETSC ERROR: or try <a class="gmail-m_-1435807948543973791gmail-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.8<wbr>.2/src/mat/impls/aij/mpi/mpipt<wbr>ap.c<br>
[0]PETSC ERROR: [0] MatView line 949
/home/lanthale/Progs/petsc-3.8<wbr>.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_-1435807948543973791gmail-m_349307368681188326moz-txt-link-freetext" href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/d<wbr>ocumentation/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>
</blockquote>
<br>
</div></div></div>
</blockquote></div><br></div></div>