Actual source code: ex48.c
2: static char help[] = "Tests various routines in MatSeqBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat A, B, C, D, Fact;
9: Vec xx, s1, s2, yy;
10: PetscInt m = 45, rows[2], cols[2], bs = 1, i, row, col, *idx, M;
11: PetscScalar rval, vals1[4], vals2[4];
12: PetscRandom rdm;
13: IS is1, is2;
14: PetscReal s1norm, s2norm, rnorm, tol = 1.e-4;
15: PetscBool flg;
16: MatFactorInfo info;
19: PetscInitialize(&argc, &args, (char *)0, help);
20: /* Test MatSetValues() and MatGetValues() */
21: PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL);
22: PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL);
23: M = m * bs;
24: MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A);
25: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
26: MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B);
27: MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
28: PetscRandomCreate(PETSC_COMM_SELF, &rdm);
29: PetscRandomSetFromOptions(rdm);
30: VecCreateSeq(PETSC_COMM_SELF, M, &xx);
31: VecDuplicate(xx, &s1);
32: VecDuplicate(xx, &s2);
33: VecDuplicate(xx, &yy);
35: /* For each row add at least 15 elements */
36: for (row = 0; row < M; row++) {
37: for (i = 0; i < 25 * bs; i++) {
38: PetscRandomGetValue(rdm, &rval);
39: col = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
40: MatSetValues(A, 1, &row, 1, &col, &rval, INSERT_VALUES);
41: MatSetValues(B, 1, &row, 1, &col, &rval, INSERT_VALUES);
42: }
43: }
45: /* Now set blocks of values */
46: for (i = 0; i < 20 * bs; i++) {
47: PetscRandomGetValue(rdm, &rval);
48: cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
49: vals1[0] = rval;
50: PetscRandomGetValue(rdm, &rval);
51: cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
52: vals1[1] = rval;
53: PetscRandomGetValue(rdm, &rval);
54: rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
55: vals1[2] = rval;
56: PetscRandomGetValue(rdm, &rval);
57: rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
58: vals1[3] = rval;
59: MatSetValues(A, 2, rows, 2, cols, vals1, INSERT_VALUES);
60: MatSetValues(B, 2, rows, 2, cols, vals1, INSERT_VALUES);
61: }
63: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
65: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
68: /* Test MatNorm() */
69: MatNorm(A, NORM_FROBENIUS, &s1norm);
70: MatNorm(B, NORM_FROBENIUS, &s2norm);
71: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
72: if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
73: MatNorm(A, NORM_INFINITY, &s1norm);
74: MatNorm(B, NORM_INFINITY, &s2norm);
75: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
76: if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
77: MatNorm(A, NORM_1, &s1norm);
78: MatNorm(B, NORM_1, &s2norm);
79: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
80: if (rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
82: /* MatShift() */
83: rval = 10 * s1norm;
84: MatShift(A, rval);
85: MatShift(B, rval);
87: /* Test MatTranspose() */
88: MatTranspose(A, MAT_INPLACE_MATRIX, &A);
89: MatTranspose(B, MAT_INPLACE_MATRIX, &B);
91: /* Now do MatGetValues() */
92: for (i = 0; i < 30; i++) {
93: PetscRandomGetValue(rdm, &rval);
94: cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
95: PetscRandomGetValue(rdm, &rval);
96: cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
97: PetscRandomGetValue(rdm, &rval);
98: rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
99: PetscRandomGetValue(rdm, &rval);
100: rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
101: MatGetValues(A, 2, rows, 2, cols, vals1);
102: MatGetValues(B, 2, rows, 2, cols, vals2);
103: PetscArraycmp(vals1, vals2, 4, &flg);
104: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatGetValues bs = %" PetscInt_FMT "\n", bs);
105: }
107: /* Test MatMult(), MatMultAdd() */
108: for (i = 0; i < 40; i++) {
109: VecSetRandom(xx, rdm);
110: VecSet(s2, 0.0);
111: MatMult(A, xx, s1);
112: MatMultAdd(A, xx, s2, s2);
113: VecNorm(s1, NORM_2, &s1norm);
114: VecNorm(s2, NORM_2, &s2norm);
115: rnorm = s2norm - s1norm;
116: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
117: }
119: /* Test MatMult() */
120: MatMultEqual(A, B, 10, &flg);
121: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMult()\n");
123: /* Test MatMultAdd() */
124: MatMultAddEqual(A, B, 10, &flg);
125: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultAdd()\n");
127: /* Test MatMultTranspose() */
128: MatMultTransposeEqual(A, B, 10, &flg);
129: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTranspose()\n");
131: /* Test MatMultTransposeAdd() */
132: MatMultTransposeAddEqual(A, B, 10, &flg);
133: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTransposeAdd()\n");
135: /* Test MatMatMult() */
136: MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &C);
137: MatSetRandom(C, rdm);
138: MatMatMult(A, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D);
139: MatMatMultEqual(A, C, D, 40, &flg);
140: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n");
141: MatDestroy(&D);
142: MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &D);
143: MatMatMult(A, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &D);
144: MatMatMultEqual(A, C, D, 40, &flg);
145: if (!flg) PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n");
147: /* Do LUFactor() on both the matrices */
148: PetscMalloc1(M, &idx);
149: for (i = 0; i < M; i++) idx[i] = i;
150: ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is1);
151: ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is2);
152: PetscFree(idx);
153: ISSetPermutation(is1);
154: ISSetPermutation(is2);
156: MatFactorInfoInitialize(&info);
158: info.fill = 2.0;
159: info.dtcol = 0.0;
160: info.zeropivot = 1.e-14;
161: info.pivotinblocks = 1.0;
163: if (bs < 4) {
164: MatGetFactor(A, "petsc", MAT_FACTOR_LU, &Fact);
165: MatLUFactorSymbolic(Fact, A, is1, is2, &info);
166: MatLUFactorNumeric(Fact, A, &info);
167: VecSetRandom(yy, rdm);
168: MatForwardSolve(Fact, yy, xx);
169: MatBackwardSolve(Fact, xx, s1);
170: MatDestroy(&Fact);
171: VecScale(s1, -1.0);
172: MatMultAdd(A, s1, yy, yy);
173: VecNorm(yy, NORM_2, &rnorm);
174: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n", (double)rnorm, bs);
175: }
177: MatLUFactor(B, is1, is2, &info);
178: MatLUFactor(A, is1, is2, &info);
180: /* Test MatSolveAdd() */
181: for (i = 0; i < 10; i++) {
182: VecSetRandom(xx, rdm);
183: VecSetRandom(yy, rdm);
184: MatSolveAdd(B, xx, yy, s2);
185: MatSolveAdd(A, xx, yy, s1);
186: VecNorm(s1, NORM_2, &s1norm);
187: VecNorm(s2, NORM_2, &s2norm);
188: rnorm = s2norm - s1norm;
189: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
190: }
192: /* Test MatSolveAdd() when x = A'b +x */
193: for (i = 0; i < 10; i++) {
194: VecSetRandom(xx, rdm);
195: VecSetRandom(s1, rdm);
196: VecCopy(s2, s1);
197: MatSolveAdd(B, xx, s2, s2);
198: MatSolveAdd(A, xx, s1, s1);
199: VecNorm(s1, NORM_2, &s1norm);
200: VecNorm(s2, NORM_2, &s2norm);
201: rnorm = s2norm - s1norm;
202: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
203: }
205: /* Test MatSolve() */
206: for (i = 0; i < 10; i++) {
207: VecSetRandom(xx, rdm);
208: MatSolve(B, xx, s2);
209: MatSolve(A, xx, s1);
210: VecNorm(s1, NORM_2, &s1norm);
211: VecNorm(s2, NORM_2, &s2norm);
212: rnorm = s2norm - s1norm;
213: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
214: }
216: /* Test MatSolveTranspose() */
217: if (bs < 8) {
218: for (i = 0; i < 10; i++) {
219: VecSetRandom(xx, rdm);
220: MatSolveTranspose(B, xx, s2);
221: MatSolveTranspose(A, xx, s1);
222: VecNorm(s1, NORM_2, &s1norm);
223: VecNorm(s2, NORM_2, &s2norm);
224: rnorm = s2norm - s1norm;
225: if (rnorm < -tol || rnorm > tol) PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs);
226: }
227: }
229: MatDestroy(&A);
230: MatDestroy(&B);
231: MatDestroy(&C);
232: MatDestroy(&D);
233: VecDestroy(&xx);
234: VecDestroy(&s1);
235: VecDestroy(&s2);
236: VecDestroy(&yy);
237: ISDestroy(&is1);
238: ISDestroy(&is2);
239: PetscRandomDestroy(&rdm);
240: PetscFinalize();
241: return 0;
242: }
244: /*TEST
246: test:
247: args: -mat_block_size {{1 2 3 4 5 6 7 8}}
249: TEST*/