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