Actual source code: ex40.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Tests the parallel case for MatIncreaseOverlap(). Input arguments are:\n\
      3:   -f <input_file> : file to load.  For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\
      4:   -nd <size>      : > 0  number of domains per processor \n\
      5:   -ov <overlap>   : >=0  amount of overlap between domains\n\n";

      7: #include <petscmat.h>

     11: int main(int argc,char **args)
     12: {
     14:   PetscInt       nd = 2,ov=1,i,start,m,n,end,lsize;
     15:   PetscMPIInt    rank;
     16:   PetscBool      flg;
     17:   Mat            A,B;
     18:   char           file[PETSC_MAX_PATH_LEN];
     19:   PetscViewer    fd;
     20:   IS             *is1,*is2;
     21:   PetscRandom    r;
     22:   PetscScalar    rand;

     24:   PetscInitialize(&argc,&args,(char*)0,help);
     25: #if defined(PETSC_USE_COMPLEX)
     26:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
     27: #else

     29:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
     30:   PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
     31:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must use -f filename to indicate a file containing a PETSc binary matrix");
     32:   PetscOptionsGetInt(NULL,"-nd",&nd,NULL);
     33:   PetscOptionsGetInt(NULL,"-ov",&ov,NULL);

     35:   /* Read matrix and RHS */
     36:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
     37:   MatCreate(PETSC_COMM_WORLD,&A);
     38:   MatSetType(A,MATMPIAIJ);
     39:   MatLoad(A,fd);
     40:   PetscViewerDestroy(&fd);

     42:   /* Read the matrix again as a sequential matrix */
     43:   PetscViewerBinaryOpen(PETSC_COMM_SELF,file,FILE_MODE_READ,&fd);
     44:   MatCreate(PETSC_COMM_SELF,&B);
     45:   MatSetType(B,MATSEQAIJ);
     46:   MatLoad(B,fd);
     47:   PetscViewerDestroy(&fd);

     49:   /* Create the IS corresponding to subdomains */
     50:   PetscMalloc(nd*sizeof(IS **),&is1);
     51:   PetscMalloc(nd*sizeof(IS **),&is2);

     53:   /* Create the random Index Sets */
     54:   MatGetSize(A,&m,&n);
     55:   PetscRandomCreate(PETSC_COMM_SELF,&r);
     56:   PetscRandomSetFromOptions(r);
     57:   for (i=0; i<nd; i++) {
     58:     PetscRandomGetValue(r,&rand);
     59:     start = (PetscInt)(rand*m);
     60:     PetscRandomGetValue(r,&rand);
     61:     end   = (PetscInt)(rand*m);
     62:     lsize =  end - start;
     63:     if (start > end) { start = end; lsize = -lsize;}
     64:     ISCreateStride(PETSC_COMM_SELF,lsize,start,1,is1+i);
     65:     ISCreateStride(PETSC_COMM_SELF,lsize,start,1,is2+i);
     66:   }
     67:   MatIncreaseOverlap(A,nd,is1,ov);
     68:   MatIncreaseOverlap(B,nd,is2,ov);

     72:   /* Now see if the serial and parallel case have the same answers */
     73:   for (i=0; i<nd; ++i) {
     74:     ISEqual(is1[i],is2[i],&flg);
     75:     PetscPrintf(PETSC_COMM_SELF,"proc:[%d], i=%D, flg =%d\n",rank,i,(int)flg);
     76:   }

     78:   /* Free allocated memory */
     79:   for (i=0; i<nd; ++i) {
     80:     ISDestroy(&is1[i]);
     81:     ISDestroy(&is2[i]);
     82:   }
     83:   PetscFree(is1);
     84:   PetscFree(is2);
     85:   PetscRandomDestroy(&r);
     86:   MatDestroy(&A);
     87:   MatDestroy(&B);

     89:   PetscFinalize();
     90: #endif
     91:   return 0;
     92: }