Actual source code: ex91.c
1: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";
3: #include <petscmat.h>
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;
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));
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));
31: PetscCall(PetscMalloc1(bs, &rows));
32: PetscCall(PetscMalloc1(bs, &cols));
33: PetscCall(PetscMalloc1(bs * bs, &vals));
34: PetscCall(PetscMalloc1(M, &idx));
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: }
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: }
66: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
67: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
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));
78: /* create a SeqSBAIJ matrix sA (= A) */
79: PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
80: PetscCall(MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA));
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: }
102: /* Test MatIncreaseOverlap() */
103: PetscCall(PetscMalloc1(nd, &is1));
104: PetscCall(PetscMalloc1(nd, &is2));
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: */
123: PetscCall(MatIncreaseOverlap(A, nd, is1, ov));
124: PetscCall(MatIncreaseOverlap(sA, nd, is2, ov));
126: for (i = 0; i < nd; ++i) {
127: PetscCall(ISSort(is1[i]));
128: PetscCall(ISSort(is2[i]));
129: }
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: }
136: PetscCall(MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA));
137: PetscCall(MatCreateSubMatrices(sA, nd, is2, is2, MAT_INITIAL_MATRIX, &submatsA));
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: }
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));
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: }
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));
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: }
204: /*TEST
206: test:
207: args: -ov 2
209: TEST*/