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