Actual source code: ex56.c
1: static char help[] = "Test the use of MatSetValuesBlocked(), MatZeroRows() for rectangular MatBAIJ matrix, test MatSetValuesBlocked() for MatSBAIJ matrix (-test_mat_sbaij).";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat A;
8: PetscInt bs = 3, m = 4, n = 6, i, j, val = 10, row[2], col[3], eval, rstart;
9: PetscMPIInt size, rank;
10: PetscScalar x[6][9], y[3][3], one = 1.0;
11: PetscBool flg, testsbaij = PETSC_FALSE;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &args, NULL, help));
15: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
16: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
18: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_mat_sbaij", &testsbaij));
20: if (testsbaij) {
21: PetscCall(MatCreateSBAIJ(PETSC_COMM_WORLD, bs, m * bs, n * bs, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 1, NULL, &A));
22: } else {
23: PetscCall(MatCreateBAIJ(PETSC_COMM_WORLD, bs, m * bs, n * bs, PETSC_DECIDE, PETSC_DECIDE, 1, NULL, 1, NULL, &A));
24: }
25: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
26: eval = 9;
28: PetscCall(PetscOptionsHasName(NULL, NULL, "-ass_extern", &flg));
29: if (flg && (size != 1)) rstart = m * ((rank + 1) % size);
30: else rstart = m * (rank);
32: row[0] = rstart + 0;
33: row[1] = rstart + 2;
34: col[0] = rstart + 0;
35: col[1] = rstart + 1;
36: col[2] = rstart + 3;
37: for (i = 0; i < 6; i++) {
38: for (j = 0; j < 9; j++) x[i][j] = (PetscScalar)val++;
39: }
41: PetscCall(MatSetValuesBlocked(A, 2, row, 3, col, &x[0][0], INSERT_VALUES));
43: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
44: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
46: /*
47: This option does not work for rectangular matrices
48: PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE));
49: */
51: PetscCall(MatSetValuesBlocked(A, 2, row, 3, col, &x[0][0], INSERT_VALUES));
53: /* Do another MatSetValues to test the case when only one local block is specified */
54: for (i = 0; i < 3; i++) {
55: for (j = 0; j < 3; j++) y[i][j] = (PetscScalar)(10 + i * eval + j);
56: }
57: PetscCall(MatSetValuesBlocked(A, 1, row, 1, col, &y[0][0], INSERT_VALUES));
58: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
59: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
61: PetscCall(PetscOptionsHasName(NULL, NULL, "-zero_rows", &flg));
62: if (flg) {
63: col[0] = rstart * bs + 0;
64: col[1] = rstart * bs + 1;
65: col[2] = rstart * bs + 2;
66: PetscCall(MatZeroRows(A, 3, col, one, 0, 0));
67: }
69: PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
71: PetscCall(MatDestroy(&A));
72: PetscCall(PetscFinalize());
73: return 0;
74: }
76: /*TEST
78: test:
79: filter: grep -v " MPI process"
81: test:
82: suffix: 4
83: nsize: 3
84: args: -ass_extern
85: filter: grep -v " MPI process"
87: test:
88: suffix: 5
89: nsize: 3
90: args: -ass_extern -zero_rows
91: filter: grep -v " MPI process"
93: TEST*/