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;
19: PetscInitialize(&argc, &args, (char *)0, help);
20: PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL);
21: PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL);
22: PetscOptionsGetInt(NULL, NULL, "-ov", &ov, NULL);
23: PetscOptionsGetInt(NULL, NULL, "-nd", &nd, NULL);
25: /* create a SeqBAIJ matrix A */
26: M = m * bs;
27: MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A);
28: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
29: PetscRandomCreate(PETSC_COMM_SELF, &rand);
30: PetscRandomSetFromOptions(rand);
32: PetscMalloc1(bs, &rows);
33: PetscMalloc1(bs, &cols);
34: PetscMalloc1(bs * bs, &vals);
35: 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: MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES);
48: }
49: /* second, add random blocks */
50: for (i = 0; i < 20 * bs; i++) {
51: PetscRandomGetValue(rand, &rval);
52: cols[0] = bs * (int)(PetscRealPart(rval) * m);
53: 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: PetscRandomGetValue(rand, &rval);
62: vals[j] = rval;
63: }
64: MatSetValues(A, bs, rows, bs, cols, vals, ADD_VALUES);
65: }
67: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
70: /* make A a symmetric matrix: A <- A^T + A */
71: MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans);
72: MatAXPY(A, one, Atrans, DIFFERENT_NONZERO_PATTERN);
73: MatDestroy(&Atrans);
74: MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans);
75: MatEqual(A, Atrans, &flg);
77: MatDestroy(&Atrans);
79: /* create a SeqSBAIJ matrix sA (= A) */
80: MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
81: MatConvert(A, MATSEQSBAIJ, MAT_INITIAL_MATRIX, &sA);
83: /* Test sA==A through MatMult() */
84: for (i = 0; i < nd; i++) {
85: MatGetSize(A, &mm, &nn);
86: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
87: VecDuplicate(xx, &s1);
88: VecDuplicate(xx, &s2);
89: for (j = 0; j < 3; j++) {
90: VecSetRandom(xx, rand);
91: MatMult(A, xx, s1);
92: MatMult(sA, xx, s2);
93: VecNorm(s1, NORM_2, &s1norm);
94: VecNorm(s2, NORM_2, &s2norm);
95: rnorm = s2norm - s1norm;
96: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm);
97: }
98: VecDestroy(&xx);
99: VecDestroy(&s1);
100: VecDestroy(&s2);
101: }
103: /* Test MatIncreaseOverlap() */
104: PetscMalloc1(nd, &is1);
105: PetscMalloc1(nd, &is2);
107: for (i = 0; i < nd; i++) {
108: PetscRandomGetValue(rand, &rval);
109: size = (int)(PetscRealPart(rval) * m);
110: for (j = 0; j < size; j++) {
111: 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: ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is1 + i);
116: ISCreateGeneral(PETSC_COMM_SELF, size * bs, idx, PETSC_COPY_VALUES, is2 + i);
117: }
118: /* for debugging */
119: /*
120: MatView(A,PETSC_VIEWER_STDOUT_SELF);
121: MatView(sA,PETSC_VIEWER_STDOUT_SELF);
122: */
124: MatIncreaseOverlap(A, nd, is1, ov);
125: MatIncreaseOverlap(sA, nd, is2, ov);
127: for (i = 0; i < nd; ++i) {
128: ISSort(is1[i]);
129: ISSort(is2[i]);
130: }
132: for (i = 0; i < nd; ++i) {
133: ISEqual(is1[i], is2[i], &flg);
135: }
137: MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA);
138: MatCreateSubMatrices(sA, nd, is2, is2, MAT_INITIAL_MATRIX, &submatsA);
140: /* Test MatMult() */
141: for (i = 0; i < nd; i++) {
142: MatGetSize(submatA[i], &mm, &nn);
143: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
144: VecDuplicate(xx, &s1);
145: VecDuplicate(xx, &s2);
146: for (j = 0; j < 3; j++) {
147: VecSetRandom(xx, rand);
148: MatMult(submatA[i], xx, s1);
149: MatMult(submatsA[i], xx, s2);
150: VecNorm(s1, NORM_2, &s1norm);
151: VecNorm(s2, NORM_2, &s2norm);
152: rnorm = s2norm - s1norm;
153: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm);
154: }
155: VecDestroy(&xx);
156: VecDestroy(&s1);
157: VecDestroy(&s2);
158: }
160: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
161: MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
162: MatCreateSubMatrices(sA, nd, is2, is2, MAT_REUSE_MATRIX, &submatsA);
164: /* Test MatMult() */
165: for (i = 0; i < nd; i++) {
166: MatGetSize(submatA[i], &mm, &nn);
167: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
168: VecDuplicate(xx, &s1);
169: VecDuplicate(xx, &s2);
170: for (j = 0; j < 3; j++) {
171: VecSetRandom(xx, rand);
172: MatMult(submatA[i], xx, s1);
173: MatMult(submatsA[i], xx, s2);
174: VecNorm(s1, NORM_2, &s1norm);
175: VecNorm(s2, NORM_2, &s2norm);
176: rnorm = s2norm - s1norm;
177: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm);
178: }
179: VecDestroy(&xx);
180: VecDestroy(&s1);
181: VecDestroy(&s2);
182: }
184: /* Free allocated memory */
185: for (i = 0; i < nd; ++i) {
186: ISDestroy(&is1[i]);
187: ISDestroy(&is2[i]);
188: }
189: MatDestroySubMatrices(nd, &submatA);
190: MatDestroySubMatrices(nd, &submatsA);
192: PetscFree(is1);
193: PetscFree(is2);
194: PetscFree(idx);
195: PetscFree(rows);
196: PetscFree(cols);
197: PetscFree(vals);
198: MatDestroy(&A);
199: MatDestroy(&sA);
200: PetscRandomDestroy(&rand);
201: PetscFinalize();
202: return 0;
203: }
205: /*TEST
207: test:
208: args: -ov 2
210: TEST*/