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;
18: PetscFunctionBeginUser;
19: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
20: /* Test MatSetValues() and MatGetValues() */
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_block_size", &bs, NULL));
22: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_size", &m, NULL));
23: M = m * bs;
24: PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF, bs, M, M, 1, NULL, &A));
25: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
26: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, M, M, 15, NULL, &B));
27: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
28: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rdm));
29: PetscCall(PetscRandomSetFromOptions(rdm));
30: PetscCall(VecCreateSeq(PETSC_COMM_SELF, M, &xx));
31: PetscCall(VecDuplicate(xx, &s1));
32: PetscCall(VecDuplicate(xx, &s2));
33: PetscCall(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: PetscCall(PetscRandomGetValue(rdm, &rval));
39: col = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
40: PetscCall(MatSetValues(A, 1, &row, 1, &col, &rval, INSERT_VALUES));
41: PetscCall(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: PetscCall(PetscRandomGetValue(rdm, &rval));
48: cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
49: vals1[0] = rval;
50: PetscCall(PetscRandomGetValue(rdm, &rval));
51: cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
52: vals1[1] = rval;
53: PetscCall(PetscRandomGetValue(rdm, &rval));
54: rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
55: vals1[2] = rval;
56: PetscCall(PetscRandomGetValue(rdm, &rval));
57: rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
58: vals1[3] = rval;
59: PetscCall(MatSetValues(A, 2, rows, 2, cols, vals1, INSERT_VALUES));
60: PetscCall(MatSetValues(B, 2, rows, 2, cols, vals1, INSERT_VALUES));
61: }
63: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
64: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
65: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
68: /* Test MatNorm() */
69: PetscCall(MatNorm(A, NORM_FROBENIUS, &s1norm));
70: PetscCall(MatNorm(B, NORM_FROBENIUS, &s2norm));
71: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
72: if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs));
73: PetscCall(MatNorm(A, NORM_INFINITY, &s1norm));
74: PetscCall(MatNorm(B, NORM_INFINITY, &s2norm));
75: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
76: if (rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs));
77: PetscCall(MatNorm(A, NORM_1, &s1norm));
78: PetscCall(MatNorm(B, NORM_1, &s2norm));
79: rnorm = PetscAbsReal(s2norm - s1norm) / s2norm;
80: if (rnorm > tol) PetscCall(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: PetscCall(MatShift(A, rval));
85: PetscCall(MatShift(B, rval));
87: /* Test MatTranspose() */
88: PetscCall(MatTranspose(A, MAT_INPLACE_MATRIX, &A));
89: PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
91: /* Now do MatGetValues() */
92: for (i = 0; i < 30; i++) {
93: PetscCall(PetscRandomGetValue(rdm, &rval));
94: cols[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
95: PetscCall(PetscRandomGetValue(rdm, &rval));
96: cols[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
97: PetscCall(PetscRandomGetValue(rdm, &rval));
98: rows[0] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
99: PetscCall(PetscRandomGetValue(rdm, &rval));
100: rows[1] = PetscMin(M - 1, (PetscInt)(PetscRealPart(rval) * M));
101: PetscCall(MatGetValues(A, 2, rows, 2, cols, vals1));
102: PetscCall(MatGetValues(B, 2, rows, 2, cols, vals2));
103: PetscCall(PetscArraycmp(vals1, vals2, 4, &flg));
104: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatGetValues bs = %" PetscInt_FMT "\n", bs));
105: }
107: /* Test MatMult(), MatMultAdd() */
108: for (i = 0; i < 40; i++) {
109: PetscCall(VecSetRandom(xx, rdm));
110: PetscCall(VecSet(s2, 0.0));
111: PetscCall(MatMult(A, xx, s1));
112: PetscCall(MatMultAdd(A, xx, s2, s2));
113: PetscCall(VecNorm(s1, NORM_2, &s1norm));
114: PetscCall(VecNorm(s2, NORM_2, &s2norm));
115: rnorm = s2norm - s1norm;
116: if (rnorm < -tol || rnorm > tol) PetscCall(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: PetscCall(MatMultEqual(A, B, 10, &flg));
121: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMult()\n"));
123: /* Test MatMultAdd() */
124: PetscCall(MatMultAddEqual(A, B, 10, &flg));
125: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultAdd()\n"));
127: /* Test MatMultTranspose() */
128: PetscCall(MatMultTransposeEqual(A, B, 10, &flg));
129: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTranspose()\n"));
131: /* Test MatMultTransposeAdd() */
132: PetscCall(MatMultTransposeAddEqual(A, B, 10, &flg));
133: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMultTransposeAdd()\n"));
135: /* Test MatMatMult() */
136: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &C));
137: PetscCall(MatSetRandom(C, rdm));
138: PetscCall(MatMatMult(A, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &D));
139: PetscCall(MatMatMultEqual(A, C, D, 40, &flg));
140: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n"));
141: PetscCall(MatDestroy(&D));
142: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, M, 40, NULL, &D));
143: PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &D));
144: PetscCall(MatMatMultEqual(A, C, D, 40, &flg));
145: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: MatMatMult()\n"));
147: /* Do LUFactor() on both the matrices */
148: PetscCall(PetscMalloc1(M, &idx));
149: for (i = 0; i < M; i++) idx[i] = i;
150: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is1));
151: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, M, idx, PETSC_COPY_VALUES, &is2));
152: PetscCall(PetscFree(idx));
153: PetscCall(ISSetPermutation(is1));
154: PetscCall(ISSetPermutation(is2));
156: PetscCall(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: PetscCall(MatGetFactor(A, "petsc", MAT_FACTOR_LU, &Fact));
165: PetscCall(MatLUFactorSymbolic(Fact, A, is1, is2, &info));
166: PetscCall(MatLUFactorNumeric(Fact, A, &info));
167: PetscCall(VecSetRandom(yy, rdm));
168: PetscCall(MatForwardSolve(Fact, yy, xx));
169: PetscCall(MatBackwardSolve(Fact, xx, s1));
170: PetscCall(MatDestroy(&Fact));
171: PetscCall(VecScale(s1, -1.0));
172: PetscCall(MatMultAdd(A, s1, yy, yy));
173: PetscCall(VecNorm(yy, NORM_2, &rnorm));
174: if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n", (double)rnorm, bs));
175: }
177: PetscCall(MatLUFactor(B, is1, is2, &info));
178: PetscCall(MatLUFactor(A, is1, is2, &info));
180: /* Test MatSolveAdd() */
181: for (i = 0; i < 10; i++) {
182: PetscCall(VecSetRandom(xx, rdm));
183: PetscCall(VecSetRandom(yy, rdm));
184: PetscCall(MatSolveAdd(B, xx, yy, s2));
185: PetscCall(MatSolveAdd(A, xx, yy, s1));
186: PetscCall(VecNorm(s1, NORM_2, &s1norm));
187: PetscCall(VecNorm(s2, NORM_2, &s2norm));
188: rnorm = s2norm - s1norm;
189: if (rnorm < -tol || rnorm > tol) PetscCall(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: PetscCall(VecSetRandom(xx, rdm));
195: PetscCall(VecSetRandom(s1, rdm));
196: PetscCall(VecCopy(s2, s1));
197: PetscCall(MatSolveAdd(B, xx, s2, s2));
198: PetscCall(MatSolveAdd(A, xx, s1, s1));
199: PetscCall(VecNorm(s1, NORM_2, &s1norm));
200: PetscCall(VecNorm(s2, NORM_2, &s2norm));
201: rnorm = s2norm - s1norm;
202: if (rnorm < -tol || rnorm > tol) PetscCall(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: PetscCall(VecSetRandom(xx, rdm));
208: PetscCall(MatSolve(B, xx, s2));
209: PetscCall(MatSolve(A, xx, s1));
210: PetscCall(VecNorm(s1, NORM_2, &s1norm));
211: PetscCall(VecNorm(s2, NORM_2, &s2norm));
212: rnorm = s2norm - s1norm;
213: if (rnorm < -tol || rnorm > tol) PetscCall(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: PetscCall(VecSetRandom(xx, rdm));
220: PetscCall(MatSolveTranspose(B, xx, s2));
221: PetscCall(MatSolveTranspose(A, xx, s1));
222: PetscCall(VecNorm(s1, NORM_2, &s1norm));
223: PetscCall(VecNorm(s2, NORM_2, &s2norm));
224: rnorm = s2norm - s1norm;
225: if (rnorm < -tol || rnorm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n", (double)s1norm, (double)s2norm, bs));
226: }
227: }
229: PetscCall(MatDestroy(&A));
230: PetscCall(MatDestroy(&B));
231: PetscCall(MatDestroy(&C));
232: PetscCall(MatDestroy(&D));
233: PetscCall(VecDestroy(&xx));
234: PetscCall(VecDestroy(&s1));
235: PetscCall(VecDestroy(&s2));
236: PetscCall(VecDestroy(&yy));
237: PetscCall(ISDestroy(&is1));
238: PetscCall(ISDestroy(&is2));
239: PetscCall(PetscRandomDestroy(&rdm));
240: PetscCall(PetscFinalize());
241: return 0;
242: }
244: /*TEST
246: test:
247: args: -mat_block_size {{1 2 3 4 5 6 7 8}}
249: TEST*/