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