<div dir="ltr"><div><div>I created a short code to illustrate my problem..<br>if you run the program with <br>-mat_view<br>-mat_view draw<br>-draw_pause -1<br><br></div>the first matrix you see has the correct structure while the one after Matconvert gives the wrong one.<br><br></div>I think it possibly has something to do with the way I use MatSetValuesBlocked after the convert?<br><div><div><br>#include <iostream><br>#include <string><br><br>#include "petscksp.h"<br>#include "petscao.h"<br>#include "petsc.h"<br>#include "HYPRE.h"<br><br>using namespace std;<br><br>void preassem(Mat & J, double * P, const int & n);<br><br>int main (int argc, char ** argv)<br>{<br>  MPI_Init(&argc, &argv);<br><br>  PetscInitialize(PETSC_NULL,<br>          PETSC_NULL,<br>          PETSC_NULL,<br>           PETSC_NULL);<br>  <br>  PetscViewer viewer;<br>  Mat J;<br>  PetscInt bs = 4;<br>  PetscInt n = 3;<br>  PetscInt bssq = bs * bs;<br>  double P[bssq];<br>  PetscInt local_size = bs * n;<br>  for (int i = 0; i < bssq; i++) {<br>    P[i] = 1.;<br>  }<br>  MatCreate(PETSC_COMM_WORLD, &J);<br><br>  MatSetBlockSize(J, bs); // So MatSetValuesBlocked can be used<br><br>  MatSetType(J, MATAIJ);<br><br>  PetscInt d_nnz[12] = {8, 8, 8, 8, 4, 4, 4, 4, 4, 4, 4, 4};<br>  <br>  <br>  MatSetSizes(J,<br>          local_size,<br>          local_size,<br>          PETSC_DECIDE,<br>          PETSC_DECIDE);<br>  <br>  MatSeqAIJSetPreallocation(J,<br>                PETSC_DEFAULT, <br>                d_nnz);<br><br>  preassem(J, P, n); // Shows the correct structure<br><br>  <br><br>  MatConvert(J, MATBAIJ, MAT_REUSE_MATRIX, & J);<br><br>  MatZeroEntries(J);     <br><br>  preassem(J, P, n); // Shows the incorrect structurex<br><br><br>  // Destroy PETSc objects.<br>  MatDestroy(& J);<br>  <br>  PetscViewerDestroy(& viewer);<br>  PetscFinalize();<br><br>  MPI_Finalize();<br><br>  return 0;<br>}<br><br>void preassem(Mat & J, double * P, const int & n)<br>{<br>  for (int GID = 0; GID < n; GID++) {<br>    MatSetValuesBlocked(J,<br>            1,<br>            & GID,<br>            1,<br>            & GID,<br>            P,<br>            ADD_VALUES);<br>  }<br>  int GID = 0;<br>  const int jGID = 2;<br>  <br>  MatSetValuesBlocked(J,<br>              1,<br>              & GID,<br>              1,<br>              & jGID,<br>              P,<br>              ADD_VALUES);<br>  MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);<br>  MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);<br><br>}<br><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Mar 17, 2015 at 12:43 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
> On Mar 17, 2015, at 11:40 AM, Chung-Kan Huang <<a href="mailto:ckhuangf@gmail.com">ckhuangf@gmail.com</a>> wrote:<br>
><br>
> Hi,<br>
><br>
> I wonder what I did wrong when using MATConvert...<br>
><br>
>   MatCreate(* comm_, &J_);<br>
><br>
>   PetscInt * d_nnz = NULL;<br>
>   PetscInt * o_nnz = NULL;<br>
>   PetscMalloc(local_size * sizeof(PetscInt), & d_nnz);<br>
>   PetscMalloc(local_size * sizeof(PetscInt), & o_nnz);<br>
>   for (int i = 0; i < local_size; i++) {<br>
>     d_nnz[i] = d_nnz_v[i];<br>
>     o_nnz[i] = o_nnz_v[i];<br>
>   }<br>
>   MatSetSizes(J_,<br>
>           local_size,<br>
>           local_size,<br>
>           PETSC_DECIDE,<br>
>           PETSC_DECIDE);<br>
><br>
>   MatSetBlockSize(J_, M_SIZE); // So MatSetValuesBlocked can be used<br>
>   MatSetType(J_, MATAIJ);<br>
><br>
>   if (comm_->Get_size() > 1) {<br>
<br>
</span>You don't need to have this if here; just call both MatMPIAIJSetPreallocation and MatSeqAIJSetPreallocation() and PETSc will use the correct one and ignore the other one automatically.<br>
<span class=""><br>
>     // MPI<br>
>     MatMPIAIJSetPreallocation(J_,<br>
>                   max_d_nz,<br>
>                   d_nnz,<br>
>                   max_o_nz,<br>
>                   o_nnz);<br>
><br>
>   } else {<br>
>     // Seq<br>
>     MatSeqAIJSetPreallocation(J_,<br>
>                   PETSC_DEFAULT,<br>
>                   d_nnz);<br>
>   }<br>
><br>
>   PetscFree(d_nnz);<br>
>   PetscFree(o_nnz);<br>
><br>
><br>
>   // Column oriented<br>
>   MatSetOption(J_, MAT_ROW_ORIENTED, PETSC_FALSE);<br>
>   // Test code to check if J_ has right structure<br>
>   preAssembeJ();<br>
><br>
</span>> <image.png><br>
<span class="">>   MatAssemblyBegin(J_, MAT_FINAL_ASSEMBLY);<br>
>   MatAssemblyEnd(J_, MAT_FINAL_ASSEMBLY);<br>
><br>
><br>
>  MatConvert(J_, MATBAIJ, MAT_REUSE_MATRIX, & J_);<br>
<br>
</span>   I cannot understand what is wrong. Is it producing the wrong matrix? Please send the entire code so we can run it ourselves.<br>
<br>
  Barry<br>
<span class=""><br>
><br>
>   MatZeroEntries(J_);<br>
> // Test code to check if J_ has right structure<br>
>   preAssembeJ();<br>
</span>> <image.png><br>
> Thanks,<br>
> Kan<br>
> Cheers<br>
><br>
<br>
</blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature"><p><strong>Cheers</strong></p></div>
</div>