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