[petsc-users] Setting up MPIBAIJ for MatMult

[Steena M]

Thank you, Barry. That took care of most of the issues. Two more things:

1) The code runs (MatMult executes) for block size =1. For block sizes > 1. For instance, block size = 3, throws:
 PETSC ERROR: Arguments are incompatible!
[0]PETSC ERROR: Cannot change block size 3 to 1!

Matrix setup:

ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filein,FILE_MODE_READ,&fd); 
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetType(A, MATMPIBAIJ);
ierr = MatSetBlockSize(A, bs); CHKERRQ(ierr);
ierr = MatLoad(A,fd);CHKERRQ(ierr); 
ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);


2) I'm not sure if my vector setup is ideal. Should I be using VecSetValues instead? 

  ierr = MatGetSize(A,&m,&n);
  
  ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
  ierr = VecSetType(x,VECMPI);
  ierr = VecSetSizes(x,PETSC_DECIDE,m);CHKERRQ(ierr);
  ierr = VecSetFromOptions(x);CHKERRQ(ierr); 

--------------------------------------------
On Wed, 11/12/14, Barry Smith <bsmith at mcs.anl.gov> wrote:

 Subject: Re: [petsc-users] Setting up MPIBAIJ for MatMult
 To: "Steena M" <stm8086 at yahoo.com>
 Cc: petsc-users at mcs.anl.gov
 Date: Wednesday, November 12, 2014, 11:44 AM
 
 
   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);
 > 
 > 
 > 
 >