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