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