<html>
  <head>
    <meta content="text/html; charset=utf-8" http-equiv="Content-Type">
  </head>
  <body bgcolor="#FFFFFF" text="#000000">
    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.<br>
    <br>
    Samuel<br>
    <br>
    <div class="moz-cite-prefix">On 12/05/2017 04:38 PM, Hong wrote:<br>
    </div>
    <blockquote
cite="mid:CAGCphBvPWTMWWjLH8S3VB=oCLWooEAVQVcvEoHPoDfpUQDHYdg@mail.gmail.com"
      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><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 moz-do-not-send="true"
                      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 moz-do-not-send="true"
                      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 moz-do-not-send="true"
                      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 moz-do-not-send="true"
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>
  </body>
</html>