1 static char help[] = "Creates a matrix, inserts some values, and tests MatCreateSubMatrices() and MatZeroEntries().\n\n"; 2 3 #include <petscmat.h> 4 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; 13 14 PetscFunctionBeginUser; 15 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 16 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_COMMON)); 17 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_COMMON)); 18 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)); 33 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; 42 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 50 PetscCall(ISDestroy(&irow)); 51 PetscCall(ISDestroy(&icol)); 52 53 /* test getting a reference on a submat */ 54 PetscCall(PetscObjectReference((PetscObject)submat)); 55 PetscCall(MatDestroySubMatrices(1, &submatrices)); 56 PetscCall(MatDestroy(&submat)); 57 } 58 59 /* Form submatrix with rows 2-4 and columns 4-8 */ 60 PetscCall(ISCreateStride(PETSC_COMM_SELF, 3, 2, 1, &irow)); 61 PetscCall(ISCreateStride(PETSC_COMM_SELF, 5, 4, 1, &icol)); 62 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices)); 63 submat = *submatrices; 64 65 /* Test reuse submatrices */ 66 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices)); 67 68 /* sviewer will cause the submatrices (one per processor) to be printed in the correct order */ 69 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices\n")); 70 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 71 PetscCall(MatView(submat, sviewer)); 72 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 73 PetscCall(PetscObjectReference((PetscObject)submat)); 74 PetscCall(MatDestroySubMatrices(1, &submatrices)); 75 PetscCall(MatDestroy(&submat)); 76 77 /* Form submatrix with rows 2-4 and all columns */ 78 PetscCall(ISDestroy(&icol)); 79 PetscCall(ISCreateStride(PETSC_COMM_SELF, 10, 0, 1, &icol)); 80 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_INITIAL_MATRIX, &submatrices)); 81 PetscCall(MatCreateSubMatrices(mat, 1, &irow, &icol, MAT_REUSE_MATRIX, &submatrices)); 82 submat = *submatrices; 83 84 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "\nSubmatrices with allcolumns\n")); 85 PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 86 PetscCall(MatView(submat, sviewer)); 87 PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &sviewer)); 88 89 /* Test MatDuplicate */ 90 PetscCall(MatDuplicate(submat, MAT_COPY_VALUES, &submat1)); 91 PetscCall(MatDestroy(&submat1)); 92 93 /* Zero the original matrix */ 94 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD, "Original zeroed matrix\n")); 95 PetscCall(MatZeroEntries(mat)); 96 PetscCall(MatView(mat, PETSC_VIEWER_STDOUT_WORLD)); 97 98 PetscCall(ISDestroy(&irow)); 99 PetscCall(ISDestroy(&icol)); 100 PetscCall(PetscObjectReference((PetscObject)submat)); 101 PetscCall(MatDestroySubMatrices(1, &submatrices)); 102 PetscCall(MatDestroy(&submat)); 103 PetscCall(MatDestroy(&mat)); 104 PetscCall(PetscFinalize()); 105 return 0; 106 } 107 108 /*TEST 109 110 test: 111 args: -mat_type aij 112 113 test: 114 suffix: 2 115 args: -mat_type dense 116 117 test: 118 suffix: 3 119 nsize: 3 120 args: -mat_type aij 121 122 test: 123 suffix: 4 124 nsize: 3 125 args: -mat_type dense 126 127 test: 128 suffix: 5 129 nsize: 3 130 args: -mat_type aij -test_all 131 132 TEST*/ 133