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, E, Bt, *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: MatSetBlockSize(B, bs);
29: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
30: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
31: PetscRandomSetFromOptions(rdm);
33: PetscMalloc1(bs, &rows);
34: PetscMalloc1(bs, &cols);
35: PetscMalloc1(bs * bs, &vals);
36: PetscMalloc1(M, &idx);
38: /* Now set blocks of values */
39: for (i = 0; i < 20 * bs; i++) {
40: PetscInt nr = 1, nc = 1;
41: PetscRandomGetValue(rdm, &rval);
42: cols[0] = bs * (int)(PetscRealPart(rval) * m);
43: PetscRandomGetValue(rdm, &rval);
44: rows[0] = bs * (int)(PetscRealPart(rval) * m);
45: for (j = 1; j < bs; j++) {
46: PetscRandomGetValue(rdm, &rval);
47: if (PetscRealPart(rval) > .5) rows[nr++] = rows[0] + j - 1;
48: }
49: for (j = 1; j < bs; j++) {
50: PetscRandomGetValue(rdm, &rval);
51: if (PetscRealPart(rval) > .5) cols[nc++] = cols[0] + j - 1;
52: }
54: for (j = 0; j < nr * nc; j++) {
55: PetscRandomGetValue(rdm, &rval);
56: vals[j] = rval;
57: }
58: MatSetValues(A, nr, rows, nc, cols, vals, ADD_VALUES);
59: MatSetValues(B, nr, rows, nc, cols, vals, ADD_VALUES);
60: }
62: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
63: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
64: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
67: /* Test MatConvert_SeqAIJ_Seq(S)BAIJ handles incompletely filled blocks */
68: MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &E);
69: MatDestroy(&E);
70: MatTranspose(B, MAT_INITIAL_MATRIX, &Bt);
71: MatAXPY(Bt, 1.0, B, DIFFERENT_NONZERO_PATTERN);
72: MatSetOption(Bt, MAT_SYMMETRIC, PETSC_TRUE);
73: MatConvert(Bt, MATSBAIJ, MAT_INITIAL_MATRIX, &E);
74: MatDestroy(&E);
75: MatDestroy(&Bt);
77: /* Test MatIncreaseOverlap() */
78: PetscMalloc1(nd, &is1);
79: PetscMalloc1(nd, &is2);
81: for (i = 0; i < nd; i++) {
82: PetscRandomGetValue(rdm, &rval);
83: lsize = (int)(PetscRealPart(rval) * m);
84: for (j = 0; j < lsize; j++) {
85: PetscRandomGetValue(rdm, &rval);
86: idx[j * bs] = bs * (int)(PetscRealPart(rval) * m);
87: for (k = 1; k < bs; k++) idx[j * bs + k] = idx[j * bs] + k;
88: }
89: ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is1 + i);
90: ISCreateGeneral(PETSC_COMM_SELF, lsize * bs, idx, PETSC_COPY_VALUES, is2 + i);
91: }
92: MatIncreaseOverlap(A, nd, is1, ov);
93: MatIncreaseOverlap(B, nd, is2, ov);
95: for (i = 0; i < nd; ++i) {
96: ISEqual(is1[i], is2[i], &flg);
98: }
100: for (i = 0; i < nd; ++i) {
101: ISSort(is1[i]);
102: ISSort(is2[i]);
103: }
105: MatCreateSubMatrices(A, nd, is1, is1, MAT_INITIAL_MATRIX, &submatA);
106: MatCreateSubMatrices(B, nd, is2, is2, MAT_INITIAL_MATRIX, &submatB);
108: /* Test MatMult() */
109: for (i = 0; i < nd; i++) {
110: MatGetSize(submatA[i], &mm, &nn);
111: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
112: VecDuplicate(xx, &s1);
113: VecDuplicate(xx, &s2);
114: for (j = 0; j < 3; j++) {
115: VecSetRandom(xx, rdm);
116: MatMult(submatA[i], xx, s1);
117: MatMult(submatB[i], xx, s2);
118: VecNorm(s1, NORM_2, &s1norm);
119: VecNorm(s2, NORM_2, &s2norm);
120: rnorm = s2norm - s1norm;
121: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm);
122: }
123: VecDestroy(&xx);
124: VecDestroy(&s1);
125: VecDestroy(&s2);
126: }
127: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
128: MatCreateSubMatrices(A, nd, is1, is1, MAT_REUSE_MATRIX, &submatA);
129: MatCreateSubMatrices(B, nd, is2, is2, MAT_REUSE_MATRIX, &submatB);
131: /* Test MatMult() */
132: for (i = 0; i < nd; i++) {
133: MatGetSize(submatA[i], &mm, &nn);
134: VecCreateSeq(PETSC_COMM_SELF, mm, &xx);
135: VecDuplicate(xx, &s1);
136: VecDuplicate(xx, &s2);
137: for (j = 0; j < 3; j++) {
138: VecSetRandom(xx, rdm);
139: MatMult(submatA[i], xx, s1);
140: MatMult(submatB[i], xx, s2);
141: VecNorm(s1, NORM_2, &s1norm);
142: VecNorm(s2, NORM_2, &s2norm);
143: rnorm = s2norm - s1norm;
144: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n", (double)s1norm, (double)s2norm);
145: }
146: VecDestroy(&xx);
147: VecDestroy(&s1);
148: VecDestroy(&s2);
149: }
151: /* Free allocated memory */
152: for (i = 0; i < nd; ++i) {
153: ISDestroy(&is1[i]);
154: ISDestroy(&is2[i]);
155: }
156: MatDestroySubMatrices(nd, &submatA);
157: MatDestroySubMatrices(nd, &submatB);
158: PetscFree(is1);
159: PetscFree(is2);
160: PetscFree(idx);
161: PetscFree(rows);
162: PetscFree(cols);
163: PetscFree(vals);
164: MatDestroy(&A);
165: MatDestroy(&B);
166: PetscRandomDestroy(&rdm);
167: PetscFinalize();
168: return 0;
169: }
171: /*TEST
173: test:
174: args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}
176: TEST*/