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