Actual source code: ex139.c
2: const char help[] = "Test MatCreateLocalRef()\n\n";
4: #include <petscmat.h>
6: static PetscErrorCode GetLocalRef(Mat A, IS isrow, IS iscol, Mat *B)
7: {
8: IS istmp;
10: PetscFunctionBegin;
11: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Extracting LocalRef with isrow:\n"));
12: PetscCall(ISOnComm(isrow, PETSC_COMM_WORLD, PETSC_COPY_VALUES, &istmp));
13: PetscCall(ISView(istmp, PETSC_VIEWER_STDOUT_WORLD));
14: PetscCall(ISDestroy(&istmp));
15: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Extracting LocalRef with iscol (only rank=0 shown):\n"));
16: PetscCall(ISOnComm(iscol, PETSC_COMM_WORLD, PETSC_COPY_VALUES, &istmp));
17: PetscCall(ISView(istmp, PETSC_VIEWER_STDOUT_WORLD));
18: PetscCall(ISDestroy(&istmp));
20: PetscCall(MatCreateLocalRef(A, isrow, iscol, B));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: int main(int argc, char *argv[])
25: {
26: MPI_Comm comm;
27: Mat J, B;
28: PetscInt i, j, k, l, rstart, rend, m, n, top_bs, row_bs, col_bs, nlocblocks, *idx, nrowblocks, ncolblocks, *ridx, *cidx, *icol, *irow;
29: PetscScalar *vals;
30: ISLocalToGlobalMapping brmap;
31: IS is0, is1;
32: PetscBool diag, blocked;
34: PetscFunctionBeginUser;
35: PetscCall(PetscInitialize(&argc, &argv, 0, help));
36: comm = PETSC_COMM_WORLD;
38: PetscOptionsBegin(comm, NULL, "LocalRef Test Options", NULL);
39: {
40: top_bs = 2;
41: row_bs = 2;
42: col_bs = 2;
43: diag = PETSC_FALSE;
44: blocked = PETSC_FALSE;
45: PetscCall(PetscOptionsInt("-top_bs", "Block size of top-level matrix", 0, top_bs, &top_bs, NULL));
46: PetscCall(PetscOptionsInt("-row_bs", "Block size of row map", 0, row_bs, &row_bs, NULL));
47: PetscCall(PetscOptionsInt("-col_bs", "Block size of col map", 0, col_bs, &col_bs, NULL));
48: PetscCall(PetscOptionsBool("-diag", "Extract a diagonal black", 0, diag, &diag, NULL));
49: PetscCall(PetscOptionsBool("-blocked", "Use block insertion", 0, blocked, &blocked, NULL));
50: }
51: PetscOptionsEnd();
53: PetscCall(MatCreate(comm, &J));
54: PetscCall(MatSetSizes(J, 6, 6, PETSC_DETERMINE, PETSC_DETERMINE));
55: PetscCall(MatSetBlockSize(J, top_bs));
56: PetscCall(MatSetFromOptions(J));
57: PetscCall(MatSeqBAIJSetPreallocation(J, top_bs, PETSC_DECIDE, 0));
58: PetscCall(MatMPIBAIJSetPreallocation(J, top_bs, PETSC_DECIDE, 0, PETSC_DECIDE, 0));
59: PetscCall(MatSetUp(J));
60: PetscCall(MatGetSize(J, &m, &n));
61: PetscCall(MatGetOwnershipRange(J, &rstart, &rend));
63: nlocblocks = (rend - rstart) / top_bs + 2;
64: PetscCall(PetscMalloc1(nlocblocks, &idx));
65: for (i = 0; i < nlocblocks; i++) idx[i] = (rstart / top_bs + i - 1 + m / top_bs) % (m / top_bs);
66: PetscCall(ISLocalToGlobalMappingCreate(comm, top_bs, nlocblocks, idx, PETSC_OWN_POINTER, &brmap));
67: PetscCall(PetscPrintf(comm, "Block ISLocalToGlobalMapping:\n"));
68: PetscCall(ISLocalToGlobalMappingView(brmap, PETSC_VIEWER_STDOUT_WORLD));
70: PetscCall(MatSetLocalToGlobalMapping(J, brmap, brmap));
71: PetscCall(ISLocalToGlobalMappingDestroy(&brmap));
73: /* Create index sets for local submatrix */
74: nrowblocks = (rend - rstart) / row_bs;
75: ncolblocks = (rend - rstart) / col_bs;
76: PetscCall(PetscMalloc2(nrowblocks, &ridx, ncolblocks, &cidx));
77: for (i = 0; i < nrowblocks; i++) ridx[i] = i + ((i > nrowblocks / 2) ^ !rstart);
78: for (i = 0; i < ncolblocks; i++) cidx[i] = i + ((i < ncolblocks / 2) ^ !!rstart);
79: PetscCall(ISCreateBlock(PETSC_COMM_SELF, row_bs, nrowblocks, ridx, PETSC_COPY_VALUES, &is0));
80: PetscCall(ISCreateBlock(PETSC_COMM_SELF, col_bs, ncolblocks, cidx, PETSC_COPY_VALUES, &is1));
81: PetscCall(PetscFree2(ridx, cidx));
83: if (diag) {
84: PetscCall(ISDestroy(&is1));
85: PetscCall(PetscObjectReference((PetscObject)is0));
86: is1 = is0;
87: ncolblocks = nrowblocks;
88: }
90: PetscCall(GetLocalRef(J, is0, is1, &B));
92: PetscCall(MatZeroEntries(J));
94: PetscCall(PetscMalloc3(row_bs, &irow, col_bs, &icol, row_bs * col_bs, &vals));
95: for (i = 0; i < nrowblocks; i++) {
96: for (j = 0; j < ncolblocks; j++) {
97: for (k = 0; k < row_bs; k++) {
98: for (l = 0; l < col_bs; l++) vals[k * col_bs + l] = i * 100000 + j * 1000 + k * 10 + l;
99: irow[k] = i * row_bs + k;
100: }
101: for (l = 0; l < col_bs; l++) icol[l] = j * col_bs + l;
102: if (blocked) {
103: PetscCall(MatSetValuesBlockedLocal(B, 1, &i, 1, &j, vals, ADD_VALUES));
104: } else {
105: PetscCall(MatSetValuesLocal(B, row_bs, irow, col_bs, icol, vals, ADD_VALUES));
106: }
107: }
108: }
109: PetscCall(PetscFree3(irow, icol, vals));
111: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
112: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
114: PetscCall(MatView(J, PETSC_VIEWER_STDOUT_WORLD));
116: PetscCall(ISDestroy(&is0));
117: PetscCall(ISDestroy(&is1));
118: PetscCall(MatDestroy(&B));
119: PetscCall(MatDestroy(&J));
120: PetscCall(PetscFinalize());
121: return 0;
122: }
124: /*TEST
126: test:
127: nsize: 2
128: filter: grep -v "type: mpi"
129: args: -blocked {{0 1}} -mat_type {{aij baij}}
131: TEST*/