Actual source code: ex12.c

  1: /*$Id: ex12.c,v 1.29 2001/08/07 03:03:07 balay Exp $*/

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

 6:  #include petscmat.h

  8: EXTERN int TestMatZeroRows_Basic(Mat,IS,const PetscScalar[] );
  9: EXTERN int TestMatZeroRows_with_no_allocation(Mat,IS,const PetscScalar[]);

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

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

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

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

 45:   TestMatZeroRows_Basic(A,is,PETSC_NULL);
 46:   TestMatZeroRows_Basic(A,is,&diag);

 48:   TestMatZeroRows_with_no_allocation(A,is,PETSC_NULL);
 49:   TestMatZeroRows_with_no_allocation(A,is,&diag);

 51:   MatDestroy(A);

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

 71:   TestMatZeroRows_Basic(A,is,PETSC_NULL);
 72:   TestMatZeroRows_Basic(A,is,&diag);

 74:   MatDestroy(A);
 75:   ISDestroy(is);
 76:   PetscFinalize();
 77:   return 0;
 78: }

 82: int TestMatZeroRows_Basic(Mat A,IS is,const PetscScalar diag[])
 83: {
 84:   Mat        B;
 85:   int        ierr;
 86:   PetscTruth keepzeroedrows;

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

 91:   PetscOptionsHasName(PETSC_NULL,"-keep_zeroed_rows",&keepzeroedrows);
 92:   if (keepzeroedrows) {
 93:     MatSetOption(B,MAT_KEEP_ZEROED_ROWS);
 94:   }

 96:   MatZeroRows(B,is,diag);
 97:   MatView(B,PETSC_VIEWER_STDOUT_WORLD);
 98:   MatDestroy(B);
 99:   return 0;
100: }

104: int TestMatZeroRows_with_no_allocation(Mat A,IS is,const PetscScalar diag[])
105: {
106:   Mat         B;
107:   int         ierr;

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

114:   MatZeroRows(B,is,diag);
115:   MatView(B,PETSC_VIEWER_STDOUT_WORLD);
116:   MatDestroy(B);
117:   return 0;
118: }