[petsc-users] Use MatConvert

[Barry Smith]

   We could do it but it would not be efficient. 

   It would be better if you create the BAIJ and then convert to AIJ, that can be efficient or to just use BAIJ.

   Barry

> On Mar 17, 2015, at 2:13 PM, Chung-Kan Huang <ckhuangf at gmail.com> wrote:
> 
> Do you have plan to make AIJ to BAIJ possible in the future?
> 
> On Tue, Mar 17, 2015 at 2:09 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>   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
> >
> 
> 
> 
> 
> -- 
> Cheers
>