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