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