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: 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;
20: PetscCall(PetscRandomCreate(comm, &rand));
21: PetscCall(VecSetRandom(a, rand));
22: PetscCall(VecSetRandom(b, rand));
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));
30: PetscCall(MatCreateDiagonal(a2, &A));
31: PetscCall(MatCreateDiagonal(b2, &B));
32: PetscCall(VecDestroy(&a2));
33: PetscCall(VecDestroy(&b2));
35: PetscCall(VecDuplicate(a, &a2));
36: PetscCall(VecDuplicate(b, &b2));
38: PetscCall(MatAXPY(A, 0.5, B, SAME_NONZERO_PATTERN));
39: PetscCall(VecAXPY(a, 0.5, b));
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");
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");
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");
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");
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");
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");
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));
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");
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");
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");
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");
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");
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");
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");
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");
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");
165: PetscCall(VecSet(a2, 0.5));
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");
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");
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");
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");
194: PetscCall(VecNorm(a2, NORM_INFINITY, &vec_norm));
195: PetscCheck(vec_norm < PETSC_SMALL, comm, PETSC_ERR_PLIB, "MatDiagonalGetInverse gave unexpected result");
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));
205: PetscCall(MatDestroy(&A));
206: PetscCall(MatDestroy(&B));
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: }
219: /*TEST
221: test:
222: suffix: 0
224: TEST*/