1 static char help[] = "Saves 4by4 block matrix.\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A;
8 PetscInt i, j;
9 PetscMPIInt size;
10 PetscViewer fd;
11 PetscScalar values[16], one = 1.0;
12 Vec x;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &args, NULL, help));
16 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
17 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, m "Can only run on one processor");
18
19 /*
20 Open binary file. Note that we use FILE_MODE_WRITE to indicate
21 writing to this file.
22 */
23 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "4by4", FILE_MODE_WRITE, &fd));
24
25 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD, 4, 12, 12, 0, 0, &A));
26
27 for (i = 0; i < 16; i++) values[i] = i;
28 for (i = 0; i < 4; i++) values[4 * i + i] += 5;
29 i = 0;
30 j = 0;
31 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
32
33 for (i = 0; i < 16; i++) values[i] = i;
34 i = 0;
35 j = 2;
36 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
37
38 for (i = 0; i < 16; i++) values[i] = i;
39 i = 1;
40 j = 0;
41 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
42
43 for (i = 0; i < 16; i++) values[i] = i;
44 for (i = 0; i < 4; i++) values[4 * i + i] += 6;
45 i = 1;
46 j = 1;
47 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
48
49 for (i = 0; i < 16; i++) values[i] = i;
50 i = 2;
51 j = 0;
52 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
53
54 for (i = 0; i < 16; i++) values[i] = i;
55 for (i = 0; i < 4; i++) values[4 * i + i] += 7;
56 i = 2;
57 j = 2;
58 PetscCall(MatSetValuesBlocked(A, 1, &i, 1, &j, values, INSERT_VALUES));
59
60 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
61 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
62 PetscCall(MatView(A, fd));
63 PetscCall(MatDestroy(&A));
64
65 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 12, &x));
66 PetscCall(VecSet(x, one));
67 PetscCall(VecView(x, fd));
68 PetscCall(VecDestroy(&x));
69
70 PetscCall(PetscViewerDestroy(&fd));
71 PetscCall(PetscFinalize());
72 return 0;
73 }
74