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(0); 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 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