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*/