[petsc-users] Use MatConvert

[Chung-Kan Huang]

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*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150317/700f90e3/attachment-0001.html>