[petsc-users] Setting up MPIBAIJ for MatMult

[Barry Smith]

  1) remove all the "preload" stuff

  2) MatLoad doesn't require or use the Mat preallocation stuff so remove all of that

   3) MatLoad determines the global matrix size from the binary file so do not set that before calling MatLoad

   4) All processes in the comm (in your case MPI_COMM_WORLD) need to call MatLoad

   5) You can call MatSetBlockSize() to set the block size before MatLoad()

   6) You do not need to call MatAssemblyBegin/End after a MatLoad() it would do nothing

> On Nov 12, 2014, at 1:32 PM, Steena M <stm8086 at yahoo.com> wrote:
> 
> 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);
> 
> 
> 
>