2: static char help[] = "Tests the use of MatZeroRows() for parallel matrices.\n\
  3: This example also tests the use of MatDuplicate() for both MPIAIJ and MPIBAIJ matrices";

 5:  #include petscmat.h

  7: EXTERN PetscErrorCode TestMatZeroRows_Basic(Mat,IS,PetscScalar);
  8: EXTERN PetscErrorCode TestMatZeroRows_with_no_allocation(Mat,IS,PetscScalar);

 12: int main(int argc,char **args)
 13: {
 14:   Mat            A;
 15:   PetscInt       i,j,m = 3,n,Ii,J,Imax;
 16:   PetscMPIInt    rank,size;
 18:   PetscScalar    v,diag=-4.0;
 19:   IS             is;

 21:   PetscInitialize(&argc,&args,(char *)0,help);
 22:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 23:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 24:   n = 2*size;

 26:   /* create A Square matrix for the five point stencil,YET AGAIN*/
 27:   MatCreate(PETSC_COMM_WORLD,&A);
 28:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 29:   MatSetFromOptions(A);
 30:   for (i=0; i<m; i++) {
 31:     for (j=2*rank; j<2*rank+2; j++) {
 32:       v = -1.0;  Ii = j + n*i;
 33:       if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 34:       if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 35:       if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 36:       if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 37:       v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 38:     }
 39:   }
 40:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 41:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 43:   /* Create AN IS required by MatZeroRows() */
 44:   Imax = n*rank; if (Imax>= n*m -m - 1) Imax = m*n - m - 1;
 45:   ISCreateStride(PETSC_COMM_SELF,m,Imax,1,&is);

 47:   TestMatZeroRows_Basic(A,is,0.0);
 48:   TestMatZeroRows_Basic(A,is,diag);

 50:   TestMatZeroRows_with_no_allocation(A,is,0.0);
 51:   TestMatZeroRows_with_no_allocation(A,is,diag);

 53:   MatDestroy(A);

 55:   /* Now Create a rectangular matrix with five point stencil (app) 
 56:    n+size is used so that this dimension is always divisible by size.
 57:    This way, we can always use bs = size for any number of procs */
 58:   MatCreate(PETSC_COMM_WORLD,&A);
 59:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*(n+size));
 60:   MatSetFromOptions(A);
 61:   for (i=0; i<m; i++) {
 62:     for (j=2*rank; j<2*rank+2; j++) {
 63:       v = -1.0;  Ii = j + n*i;
 64:       if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 65:       if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 66:       if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 67:       if (j<n+size-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 68:       v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 69:     }
 70:   }
 71:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 72:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 74:   TestMatZeroRows_Basic(A,is,0.0);
 75:   TestMatZeroRows_Basic(A,is,diag);

 77:   MatDestroy(A);
 78:   ISDestroy(is);
 79:   PetscFinalize();
 80:   return 0;
 81: }

 85: PetscErrorCode TestMatZeroRows_Basic(Mat A,IS is,PetscScalar diag)
 86: {
 87:   Mat            B;
 89:   PetscTruth     keepzeroedrows;

 91:   /* Now copy A into B, and test it with MatZeroRows() */
 92:   MatDuplicate(A,MAT_COPY_VALUES,&B);

 94:   PetscOptionsHasName(PETSC_NULL,"-keep_zeroed_rows",&keepzeroedrows);
 95:   if (keepzeroedrows) {
 96:     MatSetOption(B,MAT_KEEP_ZEROED_ROWS);
 97:   }

 99:   MatZeroRowsIS(B,is,diag);
100:   MatView(B,PETSC_VIEWER_STDOUT_WORLD);
101:   MatDestroy(B);
102:   return 0;
103: }

107: PetscErrorCode TestMatZeroRows_with_no_allocation(Mat A,IS is,PetscScalar diag)
108: {
109:   Mat            B;

112:   /* Now copy A into B, and test it with MatZeroRows() */
113:   MatDuplicate(A,MAT_COPY_VALUES,&B);
114:   /* Set this flag after assembly. This way, it affects only MatZeroRows() */
115:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR);

117:   MatZeroRowsIS(B,is,diag);
118:   MatView(B,PETSC_VIEWER_STDOUT_WORLD);
119:   MatDestroy(B);
120:   return 0;
121: }