1 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, B, E, Bt, *submatA, *submatB;
8 PetscInt bs = 1, m = 43, ov = 1, i, j, k, *rows, *cols, M, nd = 5, *idx, mm, nn, lsize;
9 PetscScalar *vals, rval;
10 IS *is1, *is2;
11 PetscRandom rdm;
12 Vec xx, s1, s2;
13 PetscReal s1norm, s2norm, rnorm, tol = PETSC_SQRT_MACHINE_EPSILON;
14 PetscBool flg;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &args, NULL, help));
18 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
22 M = m * bs;
23
24 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A));
25 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
26 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B));
27 PetscCall(MatSetBlockSize(B, bs));
28 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
29 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
30 PetscCall(PetscRandomSetFromOptions(rdm));
31
32 PetscCall(PetscMalloc1(bs, &rows));
33 PetscCall(PetscMalloc1(bs, &cols));
34 PetscCall(PetscMalloc1(bs * bs, &vals));
35 PetscCall(PetscMalloc1(M, &idx));
36
37 /* Now set blocks of values */
38 for (i = 0; i < 20 * bs; i++) {
39 PetscInt nr = 1, nc = 1;
40 PetscCall(PetscRandomGetValue(rdm, &rval));
41 cols[0] = bs * (int)(PetscRealPart(rval) * m);
42 PetscCall(PetscRandomGetValue(rdm, &rval));
43 rows[0] = bs * (int)(PetscRealPart(rval) * m);
44 for (j = 1; j < bs; j++) {
45 PetscCall(PetscRandomGetValue(rdm, &rval));
46 if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
47 }
48 for (j = 1; j < bs; j++) {
49 PetscCall(PetscRandomGetValue(rdm, &rval));
50 if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
51 }
52
53 for (j = 0; j < nr * nc; j++) {
54 PetscCall(PetscRandomGetValue(rdm, &rval));
55 vals[j] = rval;
56 }
57 PetscCall(MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES));
58 PetscCall(MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES));
59 }
60
61 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
62 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
63 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
64 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
65
66 /* Test MatConvert_SeqAIJ_Seq(S)BAIJ handles incompletely filled blocks */
67 PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &E));
68 PetscCall(MatDestroy(&E));
69 PetscCall(MatTranspose(B, MAT_INITIAL_MATRIX, &Bt));
70 PetscCall(MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN));
71 PetscCall(MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE));
72 PetscCall(MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E));
73 PetscCall(MatDestroy(&E));
74 PetscCall(MatDestroy(&Bt));
75
76 /* Test MatIncreaseOverlap() */
77 PetscCall(PetscMalloc1(nd, &is1));
78 PetscCall(PetscMalloc1(nd, &is2));
79
80 for (i = 0; i < nd; i++) {
81 PetscCall(PetscRandomGetValue(rdm, &rval));
82 lsize = (int)(PetscRealPart(rval) * m);
83 for (j = 0; j < lsize; j++) {
84 PetscCall(PetscRandomGetValue(rdm, &rval));
85 idx[j * bs] = bs * (int)(PetscRealPart(rval) * m);
86 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
87 }
88 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i));
89 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i));
90 }
91 PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
92 PetscCall(MatIncreaseOverlap(B, nd, is2, ov));
93
94 for (i = 0; i < nd; ++i) {
95 PetscCall(ISEqual(is1[i], is2[i], &flg));
96 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", flg =%d", i, (int)flg);
97 }
98
99 for (i = 0; i < nd; ++i) {
100 PetscCall(ISSort(is1[i]));
101 PetscCall(ISSort(is2[i]));
102 }
103
104 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
105 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB));
106
107 /* Test MatMult() */
108 for (i = 0; i < nd; i++) {
109 PetscCall(MatGetSize(submatA[i], &mm, &nn));
110 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
111 PetscCall(VecDuplicate(xx, &s1));
112 PetscCall(VecDuplicate(xx, &s2));
113 for (j = 0; j < 3; j++) {
114 PetscCall(VecSetRandom(xx, rdm));
115 PetscCall(MatMult(submatA[i], xx, s1));
116 PetscCall(MatMult(submatB[i], xx, s2));
117 PetscCall(VecNorm(s1, NORM_2, &s1norm));
118 PetscCall(VecNorm(s2, NORM_2, &s2norm));
119 rnorm = s2norm - s1norm;
120 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
121 }
122 PetscCall(VecDestroy(&xx));
123 PetscCall(VecDestroy(&s1));
124 PetscCall(VecDestroy(&s2));
125 }
126 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
127 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
128 PetscCall(MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB));
129
130 /* Test MatMult() */
131 for (i = 0; i < nd; i++) {
132 PetscCall(MatGetSize(submatA[i], &mm, &nn));
133 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
134 PetscCall(VecDuplicate(xx, &s1));
135 PetscCall(VecDuplicate(xx, &s2));
136 for (j = 0; j < 3; j++) {
137 PetscCall(VecSetRandom(xx, rdm));
138 PetscCall(MatMult(submatA[i], xx, s1));
139 PetscCall(MatMult(submatB[i], xx, s2));
140 PetscCall(VecNorm(s1, NORM_2, &s1norm));
141 PetscCall(VecNorm(s2, NORM_2, &s2norm));
142 rnorm = s2norm - s1norm;
143 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
144 }
145 PetscCall(VecDestroy(&xx));
146 PetscCall(VecDestroy(&s1));
147 PetscCall(VecDestroy(&s2));
148 }
149
150 /* Free allocated memory */
151 for (i = 0; i < nd; ++i) {
152 PetscCall(ISDestroy(&is1[i]));
153 PetscCall(ISDestroy(&is2[i]));
154 }
155 PetscCall(MatDestroySubMatrices(nd, &submatA));
156 PetscCall(MatDestroySubMatrices(nd, &submatB));
157 PetscCall(PetscFree(is1));
158 PetscCall(PetscFree(is2));
159 PetscCall(PetscFree(idx));
160 PetscCall(PetscFree(rows));
161 PetscCall(PetscFree(cols));
162 PetscCall(PetscFree(vals));
163 PetscCall(MatDestroy(&A));
164 PetscCall(MatDestroy(&B));
165 PetscCall(PetscRandomDestroy(&rdm));
166 PetscCall(PetscFinalize());
167 return 0;
168 }
169
170 /*TEST
171
172 test:
173 args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}
174 output_file: output/empty.out
175
176 TEST*/
177