<div dir="ltr">Do you have plan to make AIJ to BAIJ possible in the future?<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Mar 17, 2015 at 2:09 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"><br>
  I'm sorry, this is our fault. The code cannot correctly convert from AIJ to BAIJ for a block size > 1. The resulting BAIJ had a block size of 1 (which makes little sense in this case) hence the MatSetValuesBlocked() behaved in a way not expected. I have changed the code to generate an error in this case instead of being misleading.<br>
<span class="HOEnZb"><font color="#888888"><br>
<br>
   Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
> On Mar 17, 2015, at 1:34 PM, Chung-Kan Huang <<a href="mailto:ckhuangf@gmail.com">ckhuangf@gmail.com</a>> wrote:<br>
><br>
> 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>
> the first matrix you see has the correct structure while the one after Matconvert gives the wrong one.<br>
><br>
> I think it possibly has something to do with the way I use MatSetValuesBlocked after the convert?<br>
><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>
><br>
> On Tue, Mar 17, 2015 at 12:43 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><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>
> 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>
><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>
> > <image.png><br>
> >   MatAssemblyBegin(J_, MAT_FINAL_ASSEMBLY);<br>
> >   MatAssemblyEnd(J_, MAT_FINAL_ASSEMBLY);<br>
> ><br>
> ><br>
> >  MatConvert(J_, MATBAIJ, MAT_REUSE_MATRIX, & J_);<br>
><br>
>    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>
><br>
> ><br>
> >   MatZeroEntries(J_);<br>
> > // Test code to check if J_ has right structure<br>
> >   preAssembeJ();<br>
> > <image.png><br>
> > Thanks,<br>
> > Kan<br>
> > Cheers<br>
> ><br>
><br>
><br>
><br>
><br>
> --<br>
> Cheers<br>
><br>
<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br><div class="gmail_signature"><p><strong>Cheers</strong></p></div>
</div>