Actual source code: ex4.c

  1: static char help[] = "Creates a matrix, inserts some values, and tests MatCreateSubMatrices() and MatZeroEntries().\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **argv)
  6: {
  7:   Mat         mat, submat, submat1, *submatrices;
  8:   PetscInt    m = 10, n = 10, i = 4, tmp, rstart, rend;
  9:   IS          irow, icol;
 10:   PetscScalar value = 1.0;
 11:   PetscViewer sviewer;
 12:   PetscBool   allA = PETSC_FALSE;

 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 16:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON));
 17:   PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_COMMON));

 19:   PetscCall(MatCreate(PETSC_COMM_WORLD, &mat));
 20:   PetscCall(MatSetSizes(mat, PETSC_DECIDE, PETSC_DECIDE, m, n));
 21:   PetscCall(MatSetFromOptions(mat));
 22:   PetscCall(MatSetUp(mat));
 23:   PetscCall(MatGetOwnershipRange(mat, &rstart, &rend));
 24:   for (i = rstart; i < rend; i++) {
 25:     value = (PetscReal)i + 1;
 26:     tmp   = i % 5;
 27:     PetscCall(MatSetValues(mat, 1, &tmp, 1, &i, &value, INSERT_VALUES));
 28:   }
 29:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
 30:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
 31:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original matrix\n"));
 32:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));

 34:   /* Test MatCreateSubMatrix_XXX_All(), i.e., submatrix = A */
 35:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_all", &allA, NULL));
 36:   if (allA) {
 37:     PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &irow));
 38:     PetscCall(ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &icol));
 39:     PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
 40:     PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
 41:     submat = *submatrices;

 43:     /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
 44:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with all\n"));
 45:     PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "--------------------\n"));
 46:     PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 47:     PetscCall(MatView(submat, sviewer));
 48:     PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 49:     PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));

 51:     PetscCall(ISDestroy(&irow));
 52:     PetscCall(ISDestroy(&icol));

 54:     /* test getting a reference on a submat */
 55:     PetscCall(PetscObjectReference((PetscObject)submat));
 56:     PetscCall(MatDestroySubMatrices(1, &submatrices));
 57:     PetscCall(MatDestroy(&submat));
 58:   }

 60:   /* Form submatrix with rows 2-4 and columns 4-8 */
 61:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 2, 1, &irow));
 62:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 5, 4, 1, &icol));
 63:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
 64:   submat = *submatrices;

 66:   /* Test reuse submatrices */
 67:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));

 69:   /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */
 70:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices\n"));
 71:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 72:   PetscCall(MatView(submat, sviewer));
 73:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 74:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));
 75:   PetscCall(PetscObjectReference((PetscObject)submat));
 76:   PetscCall(MatDestroySubMatrices(1, &submatrices));
 77:   PetscCall(MatDestroy(&submat));

 79:   /* Form submatrix with rows 2-4 and all columns */
 80:   PetscCall(ISDestroy(&icol));
 81:   PetscCall(ISCreateStride(PETSC_COMM_SELF, 10, 0, 1, &icol));
 82:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices));
 83:   PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices));
 84:   submat = *submatrices;

 86:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with allcolumns\n"));
 87:   PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 88:   PetscCall(MatView(submat, sviewer));
 89:   PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer));
 90:   PetscCall(PetscViewerFlush(PETSC_VIEWER_STDOUT_WORLD));

 92:   /* Test MatDuplicate */
 93:   PetscCall(MatDuplicate(submat, MAT_COPY_VALUES, &submat1));
 94:   PetscCall(MatDestroy(&submat1));

 96:   /* Zero the original matrix */
 97:   PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original zeroed matrix\n"));
 98:   PetscCall(MatZeroEntries(mat));
 99:   PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD));

101:   PetscCall(ISDestroy(&irow));
102:   PetscCall(ISDestroy(&icol));
103:   PetscCall(PetscObjectReference((PetscObject)submat));
104:   PetscCall(MatDestroySubMatrices(1, &submatrices));
105:   PetscCall(MatDestroy(&submat));
106:   PetscCall(MatDestroy(&mat));
107:   PetscCall(PetscFinalize());
108:   return 0;
109: }

111: /*TEST

113:    test:
114:       args: -mat_type aij

116:    test:
117:       suffix: 2
118:       args: -mat_type dense

120:    test:
121:       suffix: 3
122:       nsize: 3
123:       args: -mat_type aij

125:    test:
126:       suffix: 4
127:       nsize: 3
128:       args: -mat_type dense

130:    test:
131:       suffix: 5
132:       nsize: 3
133:       args: -mat_type aij -test_all

135: TEST*/