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