Actual source code: ex77.c
2: static char help[] = "Tests the various sequential routines in MatSBAIJ format. Same as ex74.c except diagonal entries of the matrices are zeros.\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Vec x, y, b, s1, s2;
9: Mat A; /* linear system matrix */
10: Mat sA; /* symmetric part of the matrices */
11: PetscInt n, mbs = 16, bs = 1, nz = 3, prob = 2, i, j, col[3], row, Ii, J, n1;
12: const PetscInt *ip_ptr;
13: PetscScalar neg_one = -1.0, value[3], alpha = 0.1;
14: PetscMPIInt size;
15: IS ip, isrow, iscol;
16: PetscRandom rdm;
17: PetscBool reorder = PETSC_FALSE;
18: MatInfo minfo1, minfo2;
19: PetscReal norm1, norm2, tol = 10 * PETSC_SMALL;
22: PetscInitialize(&argc, &args, (char *)0, help);
23: MPI_Comm_size(PETSC_COMM_WORLD, &size);
25: PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
26: PetscOptionsGetInt(NULL, NULL, "-mbs", &mbs, NULL);
28: n = mbs * bs;
29: MatCreateSeqBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &A);
30: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
31: MatCreateSeqSBAIJ(PETSC_COMM_WORLD, bs, n, n, nz, NULL, &sA);
32: MatSetOption(sA, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
34: /* Test MatGetOwnershipRange() */
35: MatGetOwnershipRange(A, &Ii, &J);
36: MatGetOwnershipRange(sA, &i, &j);
37: if (i - Ii || j - J) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetOwnershipRange() in MatSBAIJ format\n");
39: /* Assemble matrix */
40: if (bs == 1) {
41: PetscOptionsGetInt(NULL, NULL, "-test_problem", &prob, NULL);
42: if (prob == 1) { /* tridiagonal matrix */
43: value[0] = -1.0;
44: value[1] = 2.0;
45: value[2] = -1.0;
46: for (i = 1; i < n - 1; i++) {
47: col[0] = i - 1;
48: col[1] = i;
49: col[2] = i + 1;
50: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
51: MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
52: }
53: i = n - 1;
54: col[0] = 0;
55: col[1] = n - 2;
56: col[2] = n - 1;
58: value[0] = 0.1;
59: value[1] = -1;
60: value[2] = 2;
61: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
62: MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
64: i = 0;
65: col[0] = 0;
66: col[1] = 1;
67: col[2] = n - 1;
69: value[0] = 2.0;
70: value[1] = -1.0;
71: value[2] = 0.1;
72: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
73: MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
74: } else if (prob == 2) { /* matrix for the five point stencil */
75: n1 = (PetscInt)(PetscSqrtReal((PetscReal)n) + 0.001);
77: for (i = 0; i < n1; i++) {
78: for (j = 0; j < n1; j++) {
79: Ii = j + n1 * i;
80: if (i > 0) {
81: J = Ii - n1;
82: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
83: MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
84: }
85: if (i < n1 - 1) {
86: J = Ii + n1;
87: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
88: MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
89: }
90: if (j > 0) {
91: J = Ii - 1;
92: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
93: MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
94: }
95: if (j < n1 - 1) {
96: J = Ii + 1;
97: MatSetValues(A, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
98: MatSetValues(sA, 1, &Ii, 1, &J, &neg_one, INSERT_VALUES);
99: }
100: }
101: }
102: }
103: } else { /* bs > 1 */
104: #if defined(DIAGB)
105: for (block = 0; block < n / bs; block++) {
106: /* diagonal blocks */
107: value[0] = -1.0;
108: value[1] = 4.0;
109: value[2] = -1.0;
110: for (i = 1 + block * bs; i < bs - 1 + block * bs; i++) {
111: col[0] = i - 1;
112: col[1] = i;
113: col[2] = i + 1;
114: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
115: MatSetValues(sA, 1, &i, 3, col, value, INSERT_VALUES);
116: }
117: i = bs - 1 + block * bs;
118: col[0] = bs - 2 + block * bs;
119: col[1] = bs - 1 + block * bs;
121: value[0] = -1.0;
122: value[1] = 4.0;
123: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
124: MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES);
126: i = 0 + block * bs;
127: col[0] = 0 + block * bs;
128: col[1] = 1 + block * bs;
130: value[0] = 4.0;
131: value[1] = -1.0;
132: MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES);
133: MatSetValues(sA, 1, &i, 2, col, value, INSERT_VALUES);
134: }
135: #endif
136: /* off-diagonal blocks */
137: value[0] = -1.0;
138: for (i = 0; i < (n / bs - 1) * bs; i++) {
139: col[0] = i + bs;
140: MatSetValues(A, 1, &i, 1, col, value, INSERT_VALUES);
141: MatSetValues(sA, 1, &i, 1, col, value, INSERT_VALUES);
142: col[0] = i;
143: row = i + bs;
144: MatSetValues(A, 1, &row, 1, col, value, INSERT_VALUES);
145: MatSetValues(sA, 1, &row, 1, col, value, INSERT_VALUES);
146: }
147: }
148: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
151: MatAssemblyBegin(sA, MAT_FINAL_ASSEMBLY);
152: MatAssemblyEnd(sA, MAT_FINAL_ASSEMBLY);
154: /* Test MatNorm() */
155: MatNorm(A, NORM_FROBENIUS, &norm1);
156: MatNorm(sA, NORM_FROBENIUS, &norm2);
157: norm1 -= norm2;
158: if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), fnorm1-fnorm2=%16.14e\n", (double)norm1);
159: MatNorm(A, NORM_INFINITY, &norm1);
160: MatNorm(sA, NORM_INFINITY, &norm2);
161: norm1 -= norm2;
162: if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm(), inf_norm1-inf_norm2=%16.14e\n", (double)norm1);
164: /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */
165: MatGetInfo(A, MAT_LOCAL, &minfo1);
166: MatGetInfo(sA, MAT_LOCAL, &minfo2);
167: i = (int)(minfo1.nz_used - minfo2.nz_used);
168: j = (int)(minfo1.nz_allocated - minfo2.nz_allocated);
169: if (i < 0 || j < 0) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetInfo()\n");
171: MatGetSize(A, &Ii, &J);
172: MatGetSize(sA, &i, &j);
173: if (i - Ii || j - J) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetSize()\n");
175: MatGetBlockSize(A, &Ii);
176: MatGetBlockSize(sA, &i);
177: if (i - Ii) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetBlockSize()\n");
179: /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */
180: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
181: PetscRandomSetFromOptions(rdm);
182: VecCreateSeq(PETSC_COMM_SELF, n, &x);
183: VecDuplicate(x, &s1);
184: VecDuplicate(x, &s2);
185: VecDuplicate(x, &y);
186: VecDuplicate(x, &b);
188: VecSetRandom(x, rdm);
190: MatDiagonalScale(A, x, x);
191: MatDiagonalScale(sA, x, x);
193: MatGetDiagonal(A, s1);
194: MatGetDiagonal(sA, s2);
195: VecNorm(s1, NORM_1, &norm1);
196: VecNorm(s2, NORM_1, &norm2);
197: norm1 -= norm2;
198: if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatGetDiagonal() \n");
200: MatScale(A, alpha);
201: MatScale(sA, alpha);
203: /* Test MatMult(), MatMultAdd() */
204: for (i = 0; i < 40; i++) {
205: VecSetRandom(x, rdm);
206: MatMult(A, x, s1);
207: MatMult(sA, x, s2);
208: VecNorm(s1, NORM_1, &norm1);
209: VecNorm(s2, NORM_1, &norm2);
210: norm1 -= norm2;
211: if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatMult(), MatDiagonalScale() or MatScale()\n");
212: }
214: for (i = 0; i < 40; i++) {
215: VecSetRandom(x, rdm);
216: VecSetRandom(y, rdm);
217: MatMultAdd(A, x, y, s1);
218: MatMultAdd(sA, x, y, s2);
219: VecNorm(s1, NORM_1, &norm1);
220: VecNorm(s2, NORM_1, &norm2);
221: norm1 -= norm2;
222: if (norm1 < -tol || norm1 > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatMultAdd(), MatDiagonalScale() or MatScale() \n");
223: }
225: /* Test MatReordering() */
226: MatGetOrdering(A, MATORDERINGNATURAL, &isrow, &iscol);
227: ip = isrow;
229: if (reorder) {
230: IS nip;
231: PetscInt *nip_ptr;
232: PetscMalloc1(mbs, &nip_ptr);
233: ISGetIndices(ip, &ip_ptr);
234: PetscArraycpy(nip_ptr, ip_ptr, mbs);
235: i = nip_ptr[1];
236: nip_ptr[1] = nip_ptr[mbs - 2];
237: nip_ptr[mbs - 2] = i;
238: i = nip_ptr[0];
239: nip_ptr[0] = nip_ptr[mbs - 1];
240: nip_ptr[mbs - 1] = i;
241: ISRestoreIndices(ip, &ip_ptr);
242: ISCreateGeneral(PETSC_COMM_SELF, mbs, nip_ptr, PETSC_COPY_VALUES, &nip);
243: PetscFree(nip_ptr);
245: MatReorderingSeqSBAIJ(sA, ip);
246: ISDestroy(&nip);
247: }
249: ISDestroy(&iscol);
250: ISDestroy(&isrow);
251: MatDestroy(&A);
252: MatDestroy(&sA);
253: VecDestroy(&x);
254: VecDestroy(&y);
255: VecDestroy(&s1);
256: VecDestroy(&s2);
257: VecDestroy(&b);
258: PetscRandomDestroy(&rdm);
260: PetscFinalize();
261: return 0;
262: }
264: /*TEST
266: test:
267: args: -bs {{1 2 3 4 5 6 7 8}}
269: TEST*/