1 static char help[] = "Test combinations of scalings, shifts and get diagonal of MATSHELL\n\n"; 2 3 #include <petscmat.h> 4 5 static PetscErrorCode myMult(Mat S,Vec x,Vec y) 6 { 7 Mat A; 8 9 PetscFunctionBegin; 10 PetscCall(MatShellGetContext(S,&A)); 11 PetscCall(MatMult(A,x,y)); 12 PetscFunctionReturn(0); 13 } 14 15 static PetscErrorCode myGetDiagonal(Mat S,Vec d) 16 { 17 Mat A; 18 19 PetscFunctionBegin; 20 PetscCall(MatShellGetContext(S,&A)); 21 PetscCall(MatGetDiagonal(A,d)); 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode shiftandscale(Mat A,Vec *D) 26 { 27 Vec ll,d,rr; 28 29 PetscFunctionBegin; 30 PetscCall(MatCreateVecs(A,&ll,&rr)); 31 PetscCall(MatCreateVecs(A,&d,NULL)); 32 PetscCall(VecSetRandom(ll,NULL)); 33 PetscCall(VecSetRandom(rr,NULL)); 34 PetscCall(VecSetRandom(d,NULL)); 35 PetscCall(MatScale(A,3.0)); 36 PetscCall(MatShift(A,-4.0)); 37 PetscCall(MatScale(A,8.0)); 38 PetscCall(MatDiagonalSet(A,d,ADD_VALUES)); 39 PetscCall(MatShift(A,9.0)); 40 PetscCall(MatScale(A,8.0)); 41 PetscCall(VecSetRandom(ll,NULL)); 42 PetscCall(VecSetRandom(rr,NULL)); 43 PetscCall(MatDiagonalScale(A,ll,rr)); 44 PetscCall(MatShift(A,2.0)); 45 PetscCall(MatScale(A,11.0)); 46 PetscCall(VecSetRandom(d,NULL)); 47 PetscCall(MatDiagonalSet(A,d,ADD_VALUES)); 48 PetscCall(VecSetRandom(ll,NULL)); 49 PetscCall(VecSetRandom(rr,NULL)); 50 PetscCall(MatDiagonalScale(A,ll,rr)); 51 PetscCall(MatShift(A,5.0)); 52 PetscCall(MatScale(A,7.0)); 53 PetscCall(MatGetDiagonal(A,d)); 54 *D = d; 55 PetscCall(VecDestroy(&ll)); 56 PetscCall(VecDestroy(&rr)); 57 PetscFunctionReturn(0); 58 } 59 60 int main(int argc,char **args) 61 { 62 Mat A,Aij,B; 63 Vec Adiag,Aijdiag; 64 PetscInt m = 3; 65 PetscReal Aijnorm,Aijdiagnorm,Bnorm,dnorm; 66 67 PetscFunctionBeginUser; 68 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 69 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 70 71 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,m,7,NULL,6,NULL,&Aij)); 72 PetscCall(MatSetRandom(Aij,NULL)); 73 PetscCall(MatSetOption(Aij,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 74 75 PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,m,Aij,&A)); 76 PetscCall(MatShellSetOperation(A,MATOP_MULT,(void (*)(void)) myMult)); 77 PetscCall(MatShellSetOperation(A,MATOP_GET_DIAGONAL,(void (*)(void)) myGetDiagonal)); 78 79 PetscCall(shiftandscale(A,&Adiag)); 80 PetscCall(MatComputeOperator(A,NULL,&B)); 81 PetscCall(shiftandscale(Aij,&Aijdiag)); 82 PetscCall(MatAXPY(Aij,-1.0,B,DIFFERENT_NONZERO_PATTERN)); 83 PetscCall(MatNorm(Aij,NORM_FROBENIUS,&Aijnorm)); 84 PetscCall(MatNorm(B,NORM_FROBENIUS,&Bnorm)); 85 PetscCheck(Aijnorm/Bnorm <= 100.0*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Altered matrices do not match, norm of difference %g",(double)(Aijnorm/Bnorm)); 86 PetscCall(VecAXPY(Aijdiag,-1.0,Adiag)); 87 PetscCall(VecNorm(Adiag,NORM_2,&dnorm)); 88 PetscCall(VecNorm(Aijdiag,NORM_2,&Aijdiagnorm)); 89 PetscCheck(Aijdiagnorm/dnorm <= 100.0*PETSC_MACHINE_EPSILON,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Altered matrices diagonals do not match, norm of difference %g",(double)(Aijdiagnorm/dnorm)); 90 PetscCall(MatDestroy(&A)); 91 PetscCall(MatDestroy(&Aij)); 92 PetscCall(VecDestroy(&Adiag)); 93 PetscCall(VecDestroy(&Aijdiag)); 94 PetscCall(MatDestroy(&B)); 95 PetscCall(PetscFinalize()); 96 return 0; 97 } 98 99 /*TEST 100 101 test: 102 nsize: {{1 2 3 4}} 103 104 TEST*/ 105