Actual source code: ex235.c

  1: static char help[] = "Test combinations of scalings, shifts and get diagonal of MATSHELL\n\n";

  3: #include <petscmat.h>

  5: static PetscErrorCode myMult(Mat S, Vec x, Vec y)
  6: {
  7:   Mat A;

  9:   PetscFunctionBegin;
 10:   PetscCall(MatShellGetContext(S, &A));
 11:   PetscCall(MatMult(A, x, y));
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

 15: static PetscErrorCode myGetDiagonal(Mat S, Vec d)
 16: {
 17:   Mat A;

 19:   PetscFunctionBegin;
 20:   PetscCall(MatShellGetContext(S, &A));
 21:   PetscCall(MatGetDiagonal(A, d));
 22:   PetscFunctionReturn(PETSC_SUCCESS);
 23: }

 25: static PetscErrorCode shiftandscale(Mat A, Vec *D)
 26: {
 27:   Vec ll, d, rr;

 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: }

 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;

 67:   PetscFunctionBeginUser;
 68:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 69:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));

 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));

 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));

 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: }

 99: /*TEST

101:     test:
102:       nsize: {{1 2 3 4}}

104: TEST*/