1 static char help[] = "Tests the various routines in MatBAIJ format.\n\
2 Input arguments are:\n\
3 -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
4
5 #include <petscmat.h>
6
main(int argc,char ** args)7 int main(int argc, char **args)
8 {
9 Mat A, B, C;
10 PetscViewer va, vb, vc;
11 Vec x, y;
12 PetscInt i, j, row, m, n, ncols1, ncols2, ct, m2, n2;
13 const PetscInt *cols1, *cols2;
14 char file[PETSC_MAX_PATH_LEN];
15 PetscBool tflg;
16 PetscScalar rval;
17 const PetscScalar *vals1, *vals2;
18 PetscReal norm1, norm2, rnorm;
19 PetscRandom r;
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscInitialize(&argc, &args, NULL, help));
23 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
24
25 /* Load the matrix as AIJ format */
26 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &va));
27 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
28 PetscCall(MatSetType(A, MATSEQAIJ));
29 PetscCall(MatLoad(A, va));
30 PetscCall(PetscViewerDestroy(&va));
31
32 /* Load the matrix as BAIJ format */
33 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &vb));
34 PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
35 PetscCall(MatSetType(B, MATSEQBAIJ));
36 PetscCall(MatLoad(B, vb));
37 PetscCall(PetscViewerDestroy(&vb));
38
39 /* Load the matrix as BAIJ format */
40 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &vc));
41 PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
42 PetscCall(MatSetType(C, MATSEQBAIJ));
43 PetscCall(MatSetFromOptions(C));
44 PetscCall(MatLoad(C, vc));
45 PetscCall(PetscViewerDestroy(&vc));
46
47 PetscCall(MatGetSize(A, &m, &n));
48 PetscCall(MatGetSize(B, &m2, &n2));
49 PetscCheck(m == m2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices are of different size. Cannot run this example");
50
51 /* Test MatEqual() */
52 PetscCall(MatEqual(B, C, &tflg));
53 PetscCheck(tflg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatEqual() failed");
54
55 /* Test MatGetDiagonal() */
56 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
57 PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &y));
58
59 PetscCall(MatGetDiagonal(A, x));
60 PetscCall(MatGetDiagonal(B, y));
61
62 PetscCall(VecEqual(x, y, &tflg));
63 PetscCheck(tflg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() failed");
64
65 /* Test MatDiagonalScale() */
66 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
67 PetscCall(PetscRandomSetFromOptions(r));
68 PetscCall(VecSetRandom(x, r));
69 PetscCall(VecSetRandom(y, r));
70
71 PetscCall(MatDiagonalScale(A, x, y));
72 PetscCall(MatDiagonalScale(B, x, y));
73 PetscCall(MatMult(A, x, y));
74 PetscCall(VecNorm(y, NORM_2, &norm1));
75 PetscCall(MatMult(B, x, y));
76 PetscCall(VecNorm(y, NORM_2, &norm2));
77 rnorm = ((norm1 - norm2) * 100) / norm1;
78 PetscCheck(rnorm >= -0.1 && rnorm <= 0.01, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatDiagonalScale() failed Norm1 %g Norm2 %g", (double)norm1, (double)norm2);
79
80 /* Test MatGetRow()/ MatRestoreRow() */
81 for (ct = 0; ct < 100; ct++) {
82 PetscCall(PetscRandomGetValue(r, &rval));
83 row = (int)(rval * m);
84 PetscCall(MatGetRow(A, row, &ncols1, &cols1, &vals1));
85 PetscCall(MatGetRow(B, row, &ncols2, &cols2, &vals2));
86
87 for (i = 0, j = 0; i < ncols1 && j < ncols2; i++) {
88 while (cols2[j] != cols1[i]) j++;
89 PetscCheck(vals1[i] == vals2[j], PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - vals incorrect.");
90 }
91 PetscCheck(i >= ncols1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetRow() failed - cols incorrect");
92
93 PetscCall(MatRestoreRow(A, row, &ncols1, &cols1, &vals1));
94 PetscCall(MatRestoreRow(B, row, &ncols2, &cols2, &vals2));
95 }
96
97 PetscCall(MatDestroy(&A));
98 PetscCall(MatDestroy(&B));
99 PetscCall(MatDestroy(&C));
100 PetscCall(VecDestroy(&x));
101 PetscCall(VecDestroy(&y));
102 PetscCall(PetscRandomDestroy(&r));
103 PetscCall(PetscFinalize());
104 return 0;
105 }
106
107 /*TEST
108
109 build:
110 requires: !complex
111
112 test:
113 args: -f ${DATAFILESPATH}/matrices/cfd.1.10 -mat_block_size 5
114 requires: !complex double datafilespath !defined(PETSC_USE_64BIT_INDICES)
115 output_file: output/empty.out
116
117 TEST*/
118