2: static char help[] = "Test MatAXPY(), and illustrate how to reduce number of mallocs used during MatSetValues() calls \n\
  3:                       Matrix C is copied from ~petsc/src/ksp/ksp/examples/tutorials/ex5.c\n\n";
  4: /*
  5:   Example: ./ex132 -mat_view ::ascii_info
  6: */

  8: #include <petscmat.h>

 12: int main(int argc,char **args)
 13: {
 14:   Mat            C,C1,C2;           /* matrix */
 15:   PetscScalar    v;
 16:   PetscInt       Ii,J,Istart,Iend;
 18:   PetscInt       i,j,m = 3,n = 2;
 19:   PetscMPIInt    size,rank;
 20:   PetscBool      mat_nonsymmetric = PETSC_FALSE;
 21:   MatInfo        info;


 24:   PetscInitialize(&argc,&args,(char*)0,help);
 25:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
 26:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 27:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 28:   n    = 2*size;

 30:   /* Set flag if we are doing a nonsymmetric problem; the default is symmetric. */
 31:   PetscOptionsGetBool(NULL,"-mat_nonsym",&mat_nonsymmetric,NULL);

 33:   MatCreate(PETSC_COMM_WORLD,&C);
 34:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 35:   MatSetFromOptions(C);
 36:   MatSeqAIJSetPreallocation(C,5,NULL);

 38:   MatGetOwnershipRange(C,&Istart,&Iend);
 39:   for (Ii=Istart; Ii<Iend; Ii++) {
 40:     v = -1.0; i = Ii/n; j = Ii - i*n;
 41:     if (i>0)   {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 42:     if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 43:     if (j>0)   {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 44:     if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 45:     v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);
 46:   }

 48:   /* Make the matrix nonsymmetric if desired */
 49:   if (mat_nonsymmetric) {
 50:     for (Ii=Istart; Ii<Iend; Ii++) {
 51:       v = -1.5; i = Ii/n;
 52:       if (i>1)   {J = Ii-n-1; MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);}
 53:     }
 54:   } else {
 55:     MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
 56:     MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);
 57:   }
 58:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 61:   /* First, create C1 = 2.0*C1 + C, C1 has less non-zeros than C */
 62:   PetscPrintf(PETSC_COMM_WORLD, "\ncreate C1 = 2.0*C1 + C, C1 has less non-zeros than C \n");
 63:   MatCreate(PETSC_COMM_WORLD,&C1);
 64:   MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 65:   MatSetFromOptions(C1);
 66:   MatSeqAIJSetPreallocation(C1,1,NULL);
 67:   for (Ii=Istart; Ii<Iend; Ii++) {
 68:     MatSetValues(C1,1,&Ii,1,&Ii,&v,ADD_VALUES);
 69:   }
 70:   MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
 71:   MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);
 72:   PetscPrintf(PETSC_COMM_WORLD," MatAXPY(C1,2.0,C,DIFFERENT_NONZERO_PATTERN)...\n");
 73:   MatAXPY(C1,2.0,C,DIFFERENT_NONZERO_PATTERN);
 74:   MatGetInfo(C1,MAT_GLOBAL_SUM,&info);
 75:   PetscPrintf(PETSC_COMM_WORLD," C1: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded);

 77:   /* Secondly, create C2 = 2.0*C2 + C, C2 has non-zero pattern of C2 + C */
 78:   PetscPrintf(PETSC_COMM_WORLD, "\ncreate C2 = 2.0*C2 + C, C2 has non-zero pattern of C2 + C \n");
 79:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&C2);
 80:   /*
 81:   MatCreate(PETSC_COMM_WORLD,&C2);
 82:   MatSetSizes(C2,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 83:   MatSetFromOptions(C2);
 84:   MatSeqAIJSetPreallocation(C2,5,NULL);
 85:   */
 86:   for (Ii=Istart; Ii<Iend; Ii++) {
 87:     v    = 1.0;
 88:     MatSetValues(C2,1,&Ii,1,&Ii,&v,ADD_VALUES);
 89:   }
 90:   MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
 91:   MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
 92:   printf(" MatAXPY(C2,2.0,C,SUBSET_NONZERO_PATTERN)...\n");
 93:   MatAXPY(C2,2.0,C,SUBSET_NONZERO_PATTERN);
 94:   MatGetInfo(C2,MAT_GLOBAL_SUM,&info);
 95:   PetscPrintf(PETSC_COMM_WORLD," C2: nz_allocated = %g; nz_used = %g; nz_unneeded = %g\n",info.nz_allocated,info.nz_used, info.nz_unneeded);

 97:   MatDestroy(&C1);
 98:   MatDestroy(&C2);
 99:   MatDestroy(&C);

101:   PetscFinalize();
102:   return 0;
103: }