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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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, NULL, 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 output_file: output/empty.out 104 105 TEST*/ 106