Actual source code: ex235.c
1: static char help[] = "Test combinations of scalings, shifts and get diagonal of MATSHELL\n\n";
3: #include <petscmat.h>
5: static PetscErrorCode myMult(Mat S, Vec x, Vec y)
6: {
7: Mat A;
9: PetscFunctionBegin;
10: PetscCall(MatShellGetContext(S, &A));
11: PetscCall(MatMult(A, x, y));
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: static PetscErrorCode myGetDiagonal(Mat S, Vec d)
16: {
17: Mat A;
19: PetscFunctionBegin;
20: PetscCall(MatShellGetContext(S, &A));
21: PetscCall(MatGetDiagonal(A, d));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: static PetscErrorCode shiftandscale(Mat A, Vec *D)
26: {
27: Vec ll, d, rr;
29: PetscFunctionBegin;
30: PetscCall(MatCreateVecs(A, &ll, &rr));
31: PetscCall(MatCreateVecs(A, &d, NULL));
32: PetscCall(VecSetRandom(ll, NULL));
33: PetscCall(VecSetRandom(rr, NULL));
34: PetscCall(VecSetRandom(d, NULL));
35: PetscCall(MatScale(A, 3.0));
36: PetscCall(MatShift(A, -4.0));
37: PetscCall(MatScale(A, 8.0));
38: PetscCall(MatDiagonalSet(A, d, ADD_VALUES));
39: PetscCall(MatShift(A, 9.0));
40: PetscCall(MatScale(A, 8.0));
41: PetscCall(VecSetRandom(ll, NULL));
42: PetscCall(VecSetRandom(rr, NULL));
43: PetscCall(MatDiagonalScale(A, ll, rr));
44: PetscCall(MatShift(A, 2.0));
45: PetscCall(MatScale(A, 11.0));
46: PetscCall(VecSetRandom(d, NULL));
47: PetscCall(MatDiagonalSet(A, d, ADD_VALUES));
48: PetscCall(VecSetRandom(ll, NULL));
49: PetscCall(VecSetRandom(rr, NULL));
50: PetscCall(MatDiagonalScale(A, ll, rr));
51: PetscCall(MatShift(A, 5.0));
52: PetscCall(MatScale(A, 7.0));
53: PetscCall(MatGetDiagonal(A, d));
54: *D = d;
55: PetscCall(VecDestroy(&ll));
56: PetscCall(VecDestroy(&rr));
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: int main(int argc, char **args)
61: {
62: Mat A, Aij, B;
63: Vec Adiag, Aijdiag;
64: PetscInt m = 3;
65: PetscReal Aijnorm, Aijdiagnorm, Bnorm, dnorm;
67: PetscFunctionBeginUser;
68: PetscCall(PetscInitialize(&argc, &args, NULL, help));
69: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
71: PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, 7, NULL, 6, NULL, &Aij));
72: PetscCall(MatSetRandom(Aij, NULL));
73: PetscCall(MatSetOption(Aij, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
75: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, Aij, &A));
76: PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))myMult));
77: PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL, (void (*)(void))myGetDiagonal));
79: PetscCall(shiftandscale(A, &Adiag));
80: PetscCall(MatComputeOperator(A, NULL, &B));
81: PetscCall(shiftandscale(Aij, &Aijdiag));
82: PetscCall(MatAXPY(Aij, -1.0, B, DIFFERENT_NONZERO_PATTERN));
83: PetscCall(MatNorm(Aij, NORM_FROBENIUS, &Aijnorm));
84: PetscCall(MatNorm(B, NORM_FROBENIUS, &Bnorm));
85: PetscCheck(Aijnorm / Bnorm <= 100.0 * PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Altered matrices do not match, norm of difference %g", (double)(Aijnorm / Bnorm));
86: PetscCall(VecAXPY(Aijdiag, -1.0, Adiag));
87: PetscCall(VecNorm(Adiag, NORM_2, &dnorm));
88: PetscCall(VecNorm(Aijdiag, NORM_2, &Aijdiagnorm));
89: PetscCheck(Aijdiagnorm / dnorm <= 100.0 * PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Altered matrices diagonals do not match, norm of difference %g", (double)(Aijdiagnorm / dnorm));
90: PetscCall(MatDestroy(&A));
91: PetscCall(MatDestroy(&Aij));
92: PetscCall(VecDestroy(&Adiag));
93: PetscCall(VecDestroy(&Aijdiag));
94: PetscCall(MatDestroy(&B));
95: PetscCall(PetscFinalize());
96: return 0;
97: }
99: /*TEST
101: test:
102: nsize: {{1 2 3 4}}
104: TEST*/