1 static char help[] = "Tests converting a SBAIJ matrix to BAIJ format with MatConvert. Modified from ex55.c\n\n";
2
3 #include <petscmat.h>
4 /* Usage: ./ex141 -mat_view */
5
main(int argc,char ** args)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];
16
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!");
22
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));
42
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;
57
58 value[0] = -1.0;
59 value[1] = 4.0;
60 PetscCall(MatSetValues(C, 1, &i, 2, col, value, INSERT_VALUES));
61
62 i = 0 + block * bs;
63 col[0] = 0 + block * bs;
64 col[1] = 1 + block * bs;
65
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 }
85
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!");
90
91 PetscCall(MatDestroy(&B));
92 PetscCall(MatDestroy(&C));
93 PetscCall(PetscFinalize());
94 return 0;
95 }
96
97 /*TEST
98
99 test:
100 output_file: output/empty.out
101
102 TEST*/
103