[petsc-users] Setting up MPIBAIJ for MatMult

[Steena M]

I am trying to read in a sparse matrix file in binary format (.dat file) to set it up as MPIBAIJ for MPIMatMult. I don't think I'm doing the file reading or the matrix and vector setup correctly. Should the file loading, matrix setup, and loading be done only on rank 0 together with the vector setup?  After looking through several examples, this is the structure of my code right now:

 Vec            x,y;  
 Mat            A;  
 PetscViewer     fd;
 int rank, global_row_size, global_col_size,bs;
 PetscBool       PetscPreLoad = PETSC_FALSE;
 char           filein[PETSC_MAX_PATH_LEN]  /*binary .dat matrix file */
 PetscScalar one = 1.0;


  /* Reading relevant block size from ENV */
  char* bs_env;
  bs_env = getenv ("BLOCK_SIZE_ENV");
  bs = (PetscInt)atoi(bs_env);

 /*Reading in matrix row/cols size from ENV. Matrices are always square: global_col_size=global_row_size*/
  
  char* m_env;
  m_env = getenv ("MATRIX_ROWS_ENV");
  global_row_size = (PetscInt)atoi(m_env);
  global_col_size = global _row_size;

  PetscInitialize(&argc,&args,(char *)0,help);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 /*Matrix reading, setup, assembly*/

  if (flg) PetscPreLoad = PETSC_TRUE;
  PetscPreLoadBegin(PetscPreLoad,"Load");
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein[PetscPreLoadIt],FILE_MODE_READ,&fd);CHKERRQ(ierr);

  
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,global_row_size,global_col_size);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetType(A, MATMPIBAIJ);
  ierr = MatMPIBAIJSetPreallocation(A,bs,5,PETSC_NULL,5,PETSC_NULL);CHKERRQ(ierr);
  ierr = MatSeqBAIJSetPreallocation(A,bs,5,PETSC_NULL);CHKERRQ(ierr);

  ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr);
  ierr = PetscLogStagePush(stage);CHKERRQ(ierr);


  if (rank==0)
  {
  ierr = MatLoad(A,fd);CHKERRQ(ierr); 
  
  }

 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 ierr = PetscLogStagePop();CHKERRQ(ierr);
 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

/*Trying two types of vector assembly*/

/* Vector setup Attempt 1*/
 
 ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetSizes(x,PETSC_DECIDE,global_row_size);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr);
  ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
  ierr = VecSetFromOptions(y);CHKERRQ(ierr);
  ierr = VecSet(x,one);CHKERRQ(ierr);
 ierr = VecSet(y,one); CHKERRQ(ierr);

/*Vector setup Attempt 2*/

ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecGetSize(x,&global_row_size);CHKERRQ(ierr);

 if (rank==0)
  {
  
  for (i = 0;  i<global_row_size; i++)
  {
  ierr = VecSetValues(x,1,&i,one,INSERT_VALUES); 
  CHKERRQ(ierr);
  }
  
  }
  ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
  ierr = VecSetFromOptions(y);CHKERRQ(ierr);
  
  ierr = VecAssemblyBegin(x); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(x); CHKERRQ(ierr);



/* SpMV*/
 ierr = MatMult(A,x,y);CHKERRQ(ierr);

  ierr = VecDestroy(&x);CHKERRQ(ierr); 
  ierr = VecDestroy(&y);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);