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 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