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 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 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, (char *)0, 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, (void (*)(void))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