[petsc-users] Use MatConvert

[Barry Smith]

  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.


   Barry

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