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 CHKERRQ(MatShellGetContext(S,&A)); 11 CHKERRQ(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 CHKERRQ(MatShellGetContext(S,&A)); 21 CHKERRQ(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 CHKERRQ(MatCreateVecs(A,&ll,&rr)); 31 CHKERRQ(MatCreateVecs(A,&d,NULL)); 32 CHKERRQ(VecSetRandom(ll,NULL)); 33 CHKERRQ(VecSetRandom(rr,NULL)); 34 CHKERRQ(VecSetRandom(d,NULL)); 35 CHKERRQ(MatScale(A,3.0)); 36 CHKERRQ(MatShift(A,-4.0)); 37 CHKERRQ(MatScale(A,8.0)); 38 CHKERRQ(MatDiagonalSet(A,d,ADD_VALUES)); 39 CHKERRQ(MatShift(A,9.0)); 40 CHKERRQ(MatScale(A,8.0)); 41 CHKERRQ(VecSetRandom(ll,NULL)); 42 CHKERRQ(VecSetRandom(rr,NULL)); 43 CHKERRQ(MatDiagonalScale(A,ll,rr)); 44 CHKERRQ(MatShift(A,2.0)); 45 CHKERRQ(MatScale(A,11.0)); 46 CHKERRQ(VecSetRandom(d,NULL)); 47 CHKERRQ(MatDiagonalSet(A,d,ADD_VALUES)); 48 CHKERRQ(VecSetRandom(ll,NULL)); 49 CHKERRQ(VecSetRandom(rr,NULL)); 50 CHKERRQ(MatDiagonalScale(A,ll,rr)); 51 CHKERRQ(MatShift(A,5.0)); 52 CHKERRQ(MatScale(A,7.0)); 53 CHKERRQ(MatGetDiagonal(A,d)); 54 *D = d; 55 CHKERRQ(VecDestroy(&ll)); 56 CHKERRQ(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 PetscErrorCode ierr; 65 PetscInt m = 3; 66 PetscReal Aijnorm,Aijdiagnorm,Bnorm,dnorm; 67 68 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 69 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 70 71 CHKERRQ(MatCreateAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,m,7,NULL,6,NULL,&Aij)); 72 CHKERRQ(MatSetRandom(Aij,NULL)); 73 CHKERRQ(MatSetOption(Aij,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE)); 74 75 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,m,Aij,&A)); 76 CHKERRQ(MatShellSetOperation(A,MATOP_MULT,(void (*)(void)) myMult)); 77 CHKERRQ(MatShellSetOperation(A,MATOP_GET_DIAGONAL,(void (*)(void)) myGetDiagonal)); 78 79 CHKERRQ(shiftandscale(A,&Adiag)); 80 CHKERRQ(MatComputeOperator(A,NULL,&B)); 81 CHKERRQ(shiftandscale(Aij,&Aijdiag)); 82 CHKERRQ(MatAXPY(Aij,-1.0,B,DIFFERENT_NONZERO_PATTERN)); 83 CHKERRQ(MatNorm(Aij,NORM_FROBENIUS,&Aijnorm)); 84 CHKERRQ(MatNorm(B,NORM_FROBENIUS,&Bnorm)); 85 PetscCheckFalse(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 CHKERRQ(VecAXPY(Aijdiag,-1.0,Adiag)); 87 CHKERRQ(VecNorm(Adiag,NORM_2,&dnorm)); 88 CHKERRQ(VecNorm(Aijdiag,NORM_2,&Aijdiagnorm)); 89 PetscCheckFalse(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 CHKERRQ(MatDestroy(&A)); 91 CHKERRQ(MatDestroy(&Aij)); 92 CHKERRQ(VecDestroy(&Adiag)); 93 CHKERRQ(VecDestroy(&Aijdiag)); 94 CHKERRQ(MatDestroy(&B)); 95 ierr = PetscFinalize(); 96 return ierr; 97 } 98 99 /*TEST 100 101 test: 102 nsize: {{1 2 3 4}} 103 104 TEST*/ 105