Actual source code: ex141.c
1: static char help[] = "Tests converting a SBAIJ matrix to BAIJ format with MatConvert. Modified from ex55.c\n\n";
3: #include <petscmat.h>
4: /* Usage: ./ex141 -mat_view */
6: int main(int argc, char **args)
7: {
8: Mat C, B;
9: PetscInt i, bs = 2, mbs, m, block, d_nz = 6, col[3];
10: PetscMPIInt size;
11: char file[PETSC_MAX_PATH_LEN];
12: PetscViewer fd;
13: PetscBool equal, loadmat;
14: PetscScalar value[4];
15: PetscInt ridx[2], cidx[2];
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, NULL, help));
19: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &loadmat));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
23: /* input matrix C */
24: if (loadmat) {
25: /* Open binary file. Load a sbaij matrix, then destroy the viewer. */
26: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
27: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
28: PetscCall(MatSetType(C, MATSEQSBAIJ));
29: PetscCall(MatLoad(C, fd));
30: PetscCall(PetscViewerDestroy(&fd));
31: } else { /* Create a sbaij mat with bs>1 */
32: mbs = 8;
33: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL));
34: m = mbs * bs;
35: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
36: PetscCall(MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, m, m));
37: PetscCall(MatSetType(C, MATSBAIJ));
38: PetscCall(MatSetFromOptions(C));
39: PetscCall(MatSeqSBAIJSetPreallocation(C, bs, d_nz, NULL));
40: PetscCall(MatSetUp(C));
41: PetscCall(MatSetOption(C, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
43: for (block = 0; block < mbs; block++) {
44: /* diagonal blocks */
45: value[0] = -1.0;
46: value[1] = 4.0;
47: value[2] = -1.0;
48: for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
49: col[0] = i - 1;
50: col[1] = i;
51: col[2] = i + 1;
52: PetscCall(MatSetValues(C, 1, &i, 3, col, value, INSERT_VALUES));
53: }
54: i = bs - 1 + block * bs;
55: col[0] = bs - 2 + block * bs;
56: col[1] = bs - 1 + block * bs;
58: value[0] = -1.0;
59: value[1] = 4.0;
60: PetscCall(MatSetValues(C, 1, &i, 2, col, value, INSERT_VALUES));
62: i = 0 + block * bs;
63: col[0] = 0 + block * bs;
64: col[1] = 1 + block * bs;
66: value[0] = 4.0;
67: value[1] = -1.0;
68: PetscCall(MatSetValues(C, 1, &i, 2, col, value, INSERT_VALUES));
69: }
70: /* off-diagonal blocks */
71: value[0] = -1.0;
72: value[1] = -0.1;
73: value[2] = 0.0;
74: value[3] = -1.0; /* row-oriented */
75: for (block = 0; block < mbs - 1; block++) {
76: for (i = 0; i < bs; i++) {
77: ridx[i] = block * bs + i;
78: cidx[i] = (block + 1) * bs + i;
79: }
80: PetscCall(MatSetValues(C, bs, ridx, bs, cidx, value, INSERT_VALUES));
81: }
82: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
83: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
84: }
86: /* convert C to BAIJ format */
87: PetscCall(MatConvert(C, MATSEQBAIJ, MAT_INITIAL_MATRIX, &B));
88: PetscCall(MatMultEqual(B, C, 10, &equal));
89: PetscCheck(equal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatConvert fails!");
91: PetscCall(MatDestroy(&B));
92: PetscCall(MatDestroy(&C));
93: PetscCall(PetscFinalize());
94: return 0;
95: }
97: /*TEST
99: test:
100: output_file: output/ex141.out
102: TEST*/