1 static char help[] = "Test combinations of scalings, shifts and get diagonal of MATSHELL\n\n";
2
3 #include <petscmat.h>
4
myMult(Mat S,Vec x,Vec y)5 static PetscErrorCode myMult(Mat S, Vec x, Vec y)
6 {
7 Mat A;
8
9 PetscFunctionBegin;
10 PetscCall(MatShellGetContext(S, &A));
11 PetscCall(MatMult(A, x, y));
12 PetscFunctionReturn(PETSC_SUCCESS);
13 }
14
myGetDiagonal(Mat S,Vec d)15 static PetscErrorCode myGetDiagonal(Mat S, Vec d)
16 {
17 Mat A;
18
19 PetscFunctionBegin;
20 PetscCall(MatShellGetContext(S, &A));
21 PetscCall(MatGetDiagonal(A, d));
22 PetscFunctionReturn(PETSC_SUCCESS);
23 }
24
shiftandscale(Mat A,Vec * D)25 static PetscErrorCode shiftandscale(Mat A, Vec *D)
26 {
27 Vec ll, d, rr;
28
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 }
59
main(int argc,char ** args)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;
66
67 PetscFunctionBeginUser;
68 PetscCall(PetscInitialize(&argc, &args, NULL, help));
69 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
70
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));
74
75 PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, Aij, &A));
76 PetscCall(MatShellSetOperation(A, MATOP_MULT, (PetscErrorCodeFn *)myMult));
77 PetscCall(MatShellSetOperation(A, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)myGetDiagonal));
78
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 }
98
99 /*TEST
100
101 test:
102 nsize: {{1 2 3 4}}
103 output_file: output/empty.out
104
105 TEST*/
106