1 const char help[] = "Test MATDIAGONAL";
2
3 #include <petsc/private/petscimpl.h>
4 #include <petscmat.h>
5
main(int argc,char ** argv)6 int main(int argc, char **argv)
7 {
8 Vec a, a2, b, b2, c, c2, A_diag, A_inv_diag;
9 Mat A, B;
10 PetscInt n = 10;
11
12 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
13
14 MPI_Comm comm = PETSC_COMM_SELF;
15 PetscCall(VecCreateSeq(comm, n, &a));
16 PetscCall(VecDuplicate(a, &b));
17 PetscCall(VecDuplicate(a, &c));
18 PetscRandom rand;
19
20 PetscCall(PetscRandomCreate(comm, &rand));
21 PetscCall(VecSetRandom(a, rand));
22 PetscCall(VecSetRandom(b, rand));
23
24 PetscCall(VecDuplicate(a, &a2));
25 PetscCall(VecCopy(a, a2));
26 PetscCall(VecDuplicate(b, &b2));
27 PetscCall(VecCopy(b, b2));
28 PetscCall(VecDuplicate(c, &c2));
29
30 PetscCall(MatCreateDiagonal(a2, &A));
31 PetscCall(MatCreateDiagonal(b2, &B));
32 PetscCall(VecDestroy(&a2));
33 PetscCall(VecDestroy(&b2));
34
35 PetscCall(VecDuplicate(a, &a2));
36 PetscCall(VecDuplicate(b, &b2));
37
38 PetscCall(MatAXPY(A, 0.5, B, SAME_NONZERO_PATTERN));
39 PetscCall(VecAXPY(a, 0.5, b));
40
41 PetscReal mat_norm, vec_norm;
42 PetscCall(VecNorm(a, NORM_2, &vec_norm));
43 PetscCall(MatNorm(A, NORM_FROBENIUS, &mat_norm));
44 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match");
45
46 // For diagonal matrix, all operator norms are the max norm of the vector
47 PetscCall(VecNorm(a, NORM_INFINITY, &vec_norm));
48 PetscCall(MatNorm(A, NORM_INFINITY, &mat_norm));
49 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match");
50 PetscCall(MatNorm(A, NORM_1, &mat_norm));
51 PetscCheck(vec_norm == mat_norm, comm, PETSC_ERR_PLIB, "Norms don't match");
52
53 PetscCall(VecPointwiseMult(c, b, a));
54 PetscCall(MatMult(A, b, c2));
55 PetscCall(VecAXPY(c2, -1.0, c));
56 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm));
57 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMult is not like VecPointwiseMultiply");
58
59 PetscCall(VecPointwiseMult(c, b, a));
60 PetscCall(VecAXPY(c, 1.0, a));
61 PetscCall(MatMultAdd(A, b, a, c2));
62 PetscCall(VecAXPY(c2, -1.0, c));
63 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm));
64 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMultAdd gave unexpected value");
65
66 PetscCall(VecSet(c, 1.2));
67 PetscCall(VecSet(c2, 1.2));
68 PetscCall(VecPointwiseMult(c, b, a));
69 PetscCall(VecAXPY(c, 1.0, c2));
70 PetscCall(MatMultAdd(A, b, c2, c2));
71 PetscCall(VecAXPY(c2, -1.0, c));
72 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm));
73 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMultAdd gave unexpected value");
74
75 PetscCall(VecPointwiseDivide(c, b, a));
76 PetscCall(MatSolve(A, b, c2));
77 PetscCall(VecAXPY(c2, -1.0, c));
78 PetscCall(VecNorm(c2, NORM_INFINITY, &vec_norm));
79 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatMult is not like VecPointwiseMultiply");
80
81 Mat A_dup;
82 PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &A_dup));
83 PetscCall(MatDestroy(&A_dup));
84 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_dup));
85 PetscCall(MatGetDiagonal(A_dup, a2));
86 PetscCall(VecAXPY(a2, -1.0, a));
87 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
88 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDuplicate with MAT_COPY_VALUES did not make a duplicate vector");
89 PetscCall(MatDestroy(&A_dup));
90
91 PetscCall(MatShift(A, 1.5));
92 PetscCall(VecShift(a, 1.5));
93 PetscCall(MatGetDiagonal(A, a2));
94 PetscCall(VecAXPY(a2, -1.0, a));
95 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
96 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatShift gave different result from VecShift");
97
98 PetscCall(MatScale(A, 0.75));
99 PetscCall(VecScale(a, 0.75));
100 PetscCall(MatGetDiagonal(A, a2));
101 PetscCall(VecAXPY(a2, -1.0, a));
102 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
103 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatScale gave different result from VecScale");
104
105 PetscCall(VecPointwiseMult(a, a, b));
106 PetscCall(MatDiagonalScale(A, b, NULL));
107 PetscCall(MatGetDiagonal(A, a2));
108 PetscCall(VecAXPY(a2, -1.0, a));
109 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
110 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalScale gave unexpected result");
111
112 PetscCall(VecPointwiseMult(a, a, b));
113 PetscCall(MatDiagonalScale(A, NULL, b));
114 PetscCall(MatGetDiagonal(A, a2));
115 PetscCall(VecAXPY(a2, -1.0, a));
116 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
117 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalScale gave unexpected result");
118
119 PetscCall(VecCopy(b, a));
120 PetscCall(MatDiagonalSet(A, b, INSERT_VALUES));
121 PetscCall(MatGetDiagonal(A, a2));
122 PetscCall(VecAXPY(a2, -1.0, a));
123 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
124 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result");
125
126 PetscCall(VecSetRandom(a, rand));
127 PetscCall(VecSetRandom(b, rand));
128 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES));
129 PetscCall(VecAXPY(a, 1.0, b));
130 PetscCall(MatDiagonalSet(A, b, ADD_VALUES));
131 PetscCall(MatGetDiagonal(A, a2));
132 PetscCall(VecAXPY(a2, -1.0, a));
133 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
134 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result");
135
136 PetscCall(VecSetRandom(a, rand));
137 PetscCall(VecSetRandom(b, rand));
138 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES));
139 PetscCall(VecPointwiseMax(a, a, b));
140 PetscCall(MatDiagonalSet(A, b, MAX_VALUES));
141 PetscCall(MatGetDiagonal(A, a2));
142 PetscCall(VecAXPY(a2, -1.0, a));
143 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
144 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result");
145
146 PetscCall(VecSetRandom(a, rand));
147 PetscCall(VecSetRandom(b, rand));
148 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES));
149 PetscCall(VecPointwiseMin(a, a, b));
150 PetscCall(MatDiagonalSet(A, b, MIN_VALUES));
151 PetscCall(MatGetDiagonal(A, a2));
152 PetscCall(VecAXPY(a2, -1.0, a));
153 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
154 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result");
155
156 PetscCall(VecSetRandom(a, rand));
157 PetscCall(VecSetRandom(b, rand));
158 PetscCall(MatDiagonalSet(A, a, INSERT_VALUES));
159 PetscCall(MatDiagonalSet(A, b, NOT_SET_VALUES));
160 PetscCall(MatGetDiagonal(A, a2));
161 PetscCall(VecAXPY(a2, -1.0, a));
162 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
163 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalSet gave unexpected result");
164
165 PetscCall(VecSet(a2, 0.5));
166
167 PetscObjectState state_pre, state_post;
168 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre));
169 PetscCall(MatDiagonalGetInverseDiagonal(A, &A_inv_diag));
170 PetscCall(MatDiagonalRestoreInverseDiagonal(A, &A_inv_diag));
171 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post));
172 PetscCheck(state_pre == state_post, comm, PETSC_ERR_PLIB, "State changed on noop");
173
174 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre));
175 PetscCall(MatDiagonalGetInverseDiagonal(A, &A_inv_diag));
176 PetscCall(VecSet(A_inv_diag, 2.0));
177 PetscCall(MatDiagonalRestoreInverseDiagonal(A, &A_inv_diag));
178 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post));
179 PetscCheck(state_pre != state_post, comm, PETSC_ERR_PLIB, "State not changed on mutation");
180
181 PetscCall(PetscObjectStateGet((PetscObject)A, &state_pre));
182 PetscCall(MatDiagonalGetDiagonal(A, &A_diag));
183 PetscCall(MatDiagonalRestoreDiagonal(A, &A_diag));
184 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post));
185 PetscCheck(state_pre == state_post, comm, PETSC_ERR_PLIB, "State changed on noop");
186
187 PetscCall(MatDiagonalGetDiagonal(A, &A_diag));
188 PetscCall(VecAXPY(a2, -1.0, A_diag));
189 PetscCall(VecSet(A_diag, 1.0));
190 PetscCall(MatDiagonalRestoreDiagonal(A, &A_diag));
191 PetscCall(PetscObjectStateGet((PetscObject)A, &state_post));
192 PetscCheck(state_pre != state_post, comm, PETSC_ERR_PLIB, "State not changed on mutation");
193
194 PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
195 PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalGetInverse gave unexpected result");
196
197 PetscCall(MatZeroEntries(A));
198 PetscCall(MatNorm(A, NORM_INFINITY, &mat_norm));
199 PetscCheck(mat_norm == 0.0, comm, PETSC_ERR_PLIB, "MatZeroEntries gave unexpected result");
200 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));
201 PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INFO));
202 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));
203 PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_SELF));
204
205 PetscCall(MatDestroy(&A));
206 PetscCall(MatDestroy(&B));
207
208 PetscCall(PetscRandomDestroy(&rand));
209 PetscCall(VecDestroy(&c2));
210 PetscCall(VecDestroy(&b2));
211 PetscCall(VecDestroy(&a2));
212 PetscCall(VecDestroy(&c));
213 PetscCall(VecDestroy(&b));
214 PetscCall(VecDestroy(&a));
215 PetscCall(PetscFinalize());
216 return 0;
217 }
218
219 /*TEST
220
221 test:
222 suffix: 0
223
224 TEST*/
225