1 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Mat A, Atrans, sA, *submatA, *submatsA;
8 PetscInt bs = 1, m = 43, ov = 1, i, j, k, *rows, *cols, M, nd = 5, *idx, mm, nn;
9 PetscMPIInt size;
10 PetscScalar *vals, rval, one = 1.0;
11 IS *is1, *is2;
12 PetscRandom rand;
13 Vec xx, s1, s2;
14 PetscReal s1norm, s2norm, rnorm, tol = 10 * PETSC_SMALL;
15 PetscBool flg;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscInitialize(&argc, &args, NULL, help));
19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
21 PetscCall(PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL));
22 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL));
23
24 /* create a SeqBAIJ matrix A */
25 M = m * bs;
26 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A));
27 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
28 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rand));
29 PetscCall(PetscRandomSetFromOptions(rand));
30
31 PetscCall(PetscMalloc1(bs, &rows));
32 PetscCall(PetscMalloc1(bs, &cols));
33 PetscCall(PetscMalloc1(bs * bs, &vals));
34 PetscCall(PetscMalloc1(M, &idx));
35
36 /* Now set blocks of random values */
37 /* first, set diagonal blocks as zero */
38 for (j = 0; j < bs * bs; j++) vals[j] = 0.0;
39 for (i = 0; i < m; i++) {
40 cols[0] = i * bs;
41 rows[0] = i * bs;
42 for (j = 1; j < bs; j++) {
43 rows[j] = rows[j - 1] + 1;
44 cols[j] = cols[j - 1] + 1;
45 }
46 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
47 }
48 /* second, add random blocks */
49 for (i = 0; i < 20 * bs; i++) {
50 PetscCall(PetscRandomGetValue(rand, &rval));
51 cols[0] = bs * (int)(PetscRealPart(rval) * m);
52 PetscCall(PetscRandomGetValue(rand, &rval));
53 rows[0] = bs * (int)(PetscRealPart(rval) * m);
54 for (j = 1; j < bs; j++) {
55 rows[j] = rows[j - 1] + 1;
56 cols[j] = cols[j - 1] + 1;
57 }
58
59 for (j = 0; j < bs * bs; j++) {
60 PetscCall(PetscRandomGetValue(rand, &rval));
61 vals[j] = rval;
62 }
63 PetscCall(MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES));
64 }
65
66 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
67 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
68
69 /* make A a symmetric matrix: A <- A^T + A */
70 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
71 PetscCall(MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN));
72 PetscCall(MatDestroy(&Atrans));
73 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
74 PetscCall(MatEqual(A, Atrans, &flg));
75 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "A+A^T is non-symmetric");
76 PetscCall(MatDestroy(&Atrans));
77
78 /* create a SeqSBAIJ matrix sA (= A) */
79 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
80 PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));
81
82 /* Test sA==A through MatMult() */
83 for (i = 0; i < nd; i++) {
84 PetscCall(MatGetSize(A, &mm, &nn));
85 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
86 PetscCall(VecDuplicate(xx, &s1));
87 PetscCall(VecDuplicate(xx, &s2));
88 for (j = 0; j < 3; j++) {
89 PetscCall(VecSetRandom(xx, rand));
90 PetscCall(MatMult(A, xx, s1));
91 PetscCall(MatMult(sA, xx, s2));
92 PetscCall(VecNorm(s1, NORM_2, &s1norm));
93 PetscCall(VecNorm(s2, NORM_2, &s2norm));
94 rnorm = s2norm - s1norm;
95 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
96 }
97 PetscCall(VecDestroy(&xx));
98 PetscCall(VecDestroy(&s1));
99 PetscCall(VecDestroy(&s2));
100 }
101
102 /* Test MatIncreaseOverlap() */
103 PetscCall(PetscMalloc1(nd, &is1));
104 PetscCall(PetscMalloc1(nd, &is2));
105
106 for (i = 0; i < nd; i++) {
107 PetscCall(PetscRandomGetValue(rand, &rval));
108 size = (int)(PetscRealPart(rval) * m);
109 for (j = 0; j < size; j++) {
110 PetscCall(PetscRandomGetValue(rand, &rval));
111 idx[j * bs] = bs * (int)(PetscRealPart(rval) * m);
112 for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
113 }
114 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is1 + i));
115 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is2 + i));
116 }
117 /* for debugging */
118 /*
119 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF));
120 PetscCall(MatView(sA,PETSC_VIEWER_STDOUT_SELF));
121 */
122
123 PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
124 PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));
125
126 for (i = 0; i < nd; ++i) {
127 PetscCall(ISSort(is1[i]));
128 PetscCall(ISSort(is2[i]));
129 }
130
131 for (i = 0; i < nd; ++i) {
132 PetscCall(ISEqual(is1[i], is2[i], &flg));
133 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "i=%" PetscInt_FMT ", is1 != is2", i);
134 }
135
136 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
137 PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_INITIAL_MATRIX, &submatsA));
138
139 /* Test MatMult() */
140 for (i = 0; i < nd; i++) {
141 PetscCall(MatGetSize(submatA[i], &mm, &nn));
142 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
143 PetscCall(VecDuplicate(xx, &s1));
144 PetscCall(VecDuplicate(xx, &s2));
145 for (j = 0; j < 3; j++) {
146 PetscCall(VecSetRandom(xx, rand));
147 PetscCall(MatMult(submatA[i], xx, s1));
148 PetscCall(MatMult(submatsA[i], xx, s2));
149 PetscCall(VecNorm(s1, NORM_2, &s1norm));
150 PetscCall(VecNorm(s2, NORM_2, &s2norm));
151 rnorm = s2norm - s1norm;
152 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
153 }
154 PetscCall(VecDestroy(&xx));
155 PetscCall(VecDestroy(&s1));
156 PetscCall(VecDestroy(&s2));
157 }
158
159 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
160 PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA));
161 PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_REUSE_MATRIX, &submatsA));
162
163 /* Test MatMult() */
164 for (i = 0; i < nd; i++) {
165 PetscCall(MatGetSize(submatA[i], &mm, &nn));
166 PetscCall(VecCreateSeq(PETSC_COMM_SELF, mm, &xx));
167 PetscCall(VecDuplicate(xx, &s1));
168 PetscCall(VecDuplicate(xx, &s2));
169 for (j = 0; j < 3; j++) {
170 PetscCall(VecSetRandom(xx, rand));
171 PetscCall(MatMult(submatA[i], xx, s1));
172 PetscCall(MatMult(submatsA[i], xx, s2));
173 PetscCall(VecNorm(s1, NORM_2, &s1norm));
174 PetscCall(VecNorm(s2, NORM_2, &s2norm));
175 rnorm = s2norm - s1norm;
176 if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm));
177 }
178 PetscCall(VecDestroy(&xx));
179 PetscCall(VecDestroy(&s1));
180 PetscCall(VecDestroy(&s2));
181 }
182
183 /* Free allocated memory */
184 for (i = 0; i < nd; ++i) {
185 PetscCall(ISDestroy(&is1[i]));
186 PetscCall(ISDestroy(&is2[i]));
187 }
188 PetscCall(MatDestroySubMatrices(nd, &submatA));
189 PetscCall(MatDestroySubMatrices(nd, &submatsA));
190
191 PetscCall(PetscFree(is1));
192 PetscCall(PetscFree(is2));
193 PetscCall(PetscFree(idx));
194 PetscCall(PetscFree(rows));
195 PetscCall(PetscFree(cols));
196 PetscCall(PetscFree(vals));
197 PetscCall(MatDestroy(&A));
198 PetscCall(MatDestroy(&sA));
199 PetscCall(PetscRandomDestroy(&rand));
200 PetscCall(PetscFinalize());
201 return 0;
202 }
203
204 /*TEST
205
206 test:
207 args: -ov 2
208 output_file: output/empty.out
209
210 TEST*/
211