[petsc-users] Use MatConvert

[Chung-Kan Huang]

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