1 static char help[] = "Tests incorrect use of MatDiagonalSet() for SHELL matrices\n\n";
2
3 #include <petscmat.h>
4
5 typedef struct _n_User *User;
6 struct _n_User {
7 Mat B;
8 };
9
MatGetDiagonal_User(Mat A,Vec X)10 static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
11 {
12 User user;
13
14 PetscFunctionBegin;
15 PetscCall(MatShellGetContext(A, &user));
16 PetscCall(MatGetDiagonal(user->B, X));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
main(int argc,char ** args)20 int main(int argc, char **args)
21 {
22 const PetscScalar xvals[] = {11, 13};
23 const PetscInt inds[] = {0, 1};
24 PetscScalar avals[] = {2, 3, 5, 7};
25 Mat A, S;
26 Vec X, Y;
27 User user;
28
29 PetscFunctionBeginUser;
30 PetscCall(PetscInitialize(&argc, &args, NULL, help));
31 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A));
32 PetscCall(MatSetUp(A));
33 PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));
34 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
35 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
36 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
37 PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
38 PetscCall(VecAssemblyBegin(X));
39 PetscCall(VecAssemblyEnd(X));
40 PetscCall(VecDuplicate(X, &Y));
41
42 PetscCall(PetscNew(&user));
43 user->B = A;
44
45 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
46 PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_User));
47 PetscCall(MatSetUp(S));
48
49 PetscCall(MatShift(S, 42));
50 PetscCall(MatGetDiagonal(S, Y));
51 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
52 PetscCall(MatDiagonalSet(S, X, ADD_VALUES));
53 PetscCall(MatGetDiagonal(S, Y));
54 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
55 PetscCall(MatScale(S, 42));
56 PetscCall(MatGetDiagonal(S, Y));
57 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
58
59 PetscCall(MatDestroy(&A));
60 PetscCall(MatDestroy(&S));
61 PetscCall(VecDestroy(&X));
62 PetscCall(VecDestroy(&Y));
63 PetscCall(PetscFree(user));
64 PetscCall(PetscFinalize());
65 return 0;
66 }
67
68 /*TEST
69
70 test:
71 args: -malloc_dump
72
73 TEST*/
74