1 static char help[] = "Tests MatCopy() for SHELL matrices\n\n";
2
3 #include <petscmat.h>
4
5 typedef struct _n_User *User;
6 struct _n_User {
7 Mat A;
8 };
9
MatMult_User(Mat A,Vec X,Vec Y)10 static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
11 {
12 User user;
13
14 PetscFunctionBegin;
15 PetscCall(MatShellGetContext(A, &user));
16 PetscCall(MatMult(user->A, X, Y));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
MatCopy_User(Mat A,Mat B,MatStructure str)20 static PetscErrorCode MatCopy_User(Mat A, Mat B, MatStructure str)
21 {
22 User userA, userB;
23
24 PetscFunctionBegin;
25 PetscCall(MatShellGetContext(A, &userA));
26 if (userA) {
27 PetscCall(PetscNew(&userB));
28 PetscCall(MatDuplicate(userA->A, MAT_COPY_VALUES, &userB->A));
29 PetscCall(MatShellSetContext(B, userB));
30 }
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
MatDestroy_User(Mat A)34 static PetscErrorCode MatDestroy_User(Mat A)
35 {
36 User user;
37
38 PetscFunctionBegin;
39 PetscCall(MatShellGetContext(A, &user));
40 if (user) {
41 PetscCall(MatDestroy(&user->A));
42 PetscCall(PetscFree(user));
43 }
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
main(int argc,char ** args)47 int main(int argc, char **args)
48 {
49 const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19};
50 const PetscInt inds[] = {0, 1};
51 PetscScalar avals[] = {2, 3, 5, 7};
52 Mat S1, S2;
53 Vec X, Y;
54 User user;
55
56 PetscFunctionBeginUser;
57 PetscCall(PetscInitialize(&argc, &args, NULL, help));
58
59 PetscCall(PetscNew(&user));
60 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &user->A));
61 PetscCall(MatSetUp(user->A));
62 PetscCall(MatSetValues(user->A, 2, inds, 2, inds, avals, INSERT_VALUES));
63 PetscCall(MatAssemblyBegin(user->A, MAT_FINAL_ASSEMBLY));
64 PetscCall(MatAssemblyEnd(user->A, MAT_FINAL_ASSEMBLY));
65 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
66 PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
67 PetscCall(VecAssemblyBegin(X));
68 PetscCall(VecAssemblyEnd(X));
69 PetscCall(VecDuplicate(X, &Y));
70 PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES));
71 PetscCall(VecAssemblyBegin(Y));
72 PetscCall(VecAssemblyEnd(Y));
73
74 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S1));
75 PetscCall(MatSetUp(S1));
76 PetscCall(MatShellSetOperation(S1, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User));
77 PetscCall(MatShellSetOperation(S1, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_User));
78 PetscCall(MatShellSetOperation(S1, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_User));
79 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, NULL, &S2));
80 PetscCall(MatSetUp(S2));
81 PetscCall(MatShellSetOperation(S2, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User));
82 PetscCall(MatShellSetOperation(S2, MATOP_COPY, (PetscErrorCodeFn *)MatCopy_User));
83 PetscCall(MatShellSetOperation(S2, MATOP_DESTROY, (PetscErrorCodeFn *)MatDestroy_User));
84
85 PetscCall(MatScale(S1, 31));
86 PetscCall(MatShift(S1, 37));
87 PetscCall(MatDiagonalScale(S1, X, Y));
88 PetscCall(MatCopy(S1, S2, SAME_NONZERO_PATTERN));
89 PetscCall(MatMult(S1, X, Y));
90 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
91 PetscCall(MatMult(S2, X, Y));
92 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
93
94 PetscCall(MatDestroy(&S1));
95 PetscCall(MatDestroy(&S2));
96 PetscCall(VecDestroy(&X));
97 PetscCall(VecDestroy(&Y));
98 PetscCall(PetscFinalize());
99 return 0;
100 }
101
102 /*TEST
103
104 test:
105 args: -malloc_dump
106
107 TEST*/
108