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:   MatShellGetContext(S, &A);
 10:   MatMult(A, x, y);
 11:   return 0;
 12: }

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

 18:   MatShellGetContext(S, &A);
 19:   MatGetDiagonal(A, d);
 20:   return 0;
 21: }

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

 27:   MatCreateVecs(A, &ll, &rr);
 28:   MatCreateVecs(A, &d, NULL);
 29:   VecSetRandom(ll, NULL);
 30:   VecSetRandom(rr, NULL);
 31:   VecSetRandom(d, NULL);
 32:   MatScale(A, 3.0);
 33:   MatShift(A, -4.0);
 34:   MatScale(A, 8.0);
 35:   MatDiagonalSet(A, d, ADD_VALUES);
 36:   MatShift(A, 9.0);
 37:   MatScale(A, 8.0);
 38:   VecSetRandom(ll, NULL);
 39:   VecSetRandom(rr, NULL);
 40:   MatDiagonalScale(A, ll, rr);
 41:   MatShift(A, 2.0);
 42:   MatScale(A, 11.0);
 43:   VecSetRandom(d, NULL);
 44:   MatDiagonalSet(A, d, ADD_VALUES);
 45:   VecSetRandom(ll, NULL);
 46:   VecSetRandom(rr, NULL);
 47:   MatDiagonalScale(A, ll, rr);
 48:   MatShift(A, 5.0);
 49:   MatScale(A, 7.0);
 50:   MatGetDiagonal(A, d);
 51:   *D = d;
 52:   VecDestroy(&ll);
 53:   VecDestroy(&rr);
 54:   return 0;
 55: }

 57: int main(int argc, char **args)
 58: {
 59:   Mat       A, Aij, B;
 60:   Vec       Adiag, Aijdiag;
 61:   PetscInt  m = 3;
 62:   PetscReal Aijnorm, Aijdiagnorm, Bnorm, dnorm;

 65:   PetscInitialize(&argc, &args, (char *)0, help);
 66:   PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);

 68:   MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, 7, NULL, 6, NULL, &Aij);
 69:   MatSetRandom(Aij, NULL);
 70:   MatSetOption(Aij, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE);

 72:   MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, m, Aij, &A);
 73:   MatShellSetOperation(A, MATOP_MULT, (void (*)(void))myMult);
 74:   MatShellSetOperation(A, MATOP_GET_DIAGONAL, (void (*)(void))myGetDiagonal);

 76:   shiftandscale(A, &Adiag);
 77:   MatComputeOperator(A, NULL, &B);
 78:   shiftandscale(Aij, &Aijdiag);
 79:   MatAXPY(Aij, -1.0, B, DIFFERENT_NONZERO_PATTERN);
 80:   MatNorm(Aij, NORM_FROBENIUS, &Aijnorm);
 81:   MatNorm(B, NORM_FROBENIUS, &Bnorm);
 83:   VecAXPY(Aijdiag, -1.0, Adiag);
 84:   VecNorm(Adiag, NORM_2, &dnorm);
 85:   VecNorm(Aijdiag, NORM_2, &Aijdiagnorm);
 87:   MatDestroy(&A);
 88:   MatDestroy(&Aij);
 89:   VecDestroy(&Adiag);
 90:   VecDestroy(&Aijdiag);
 91:   MatDestroy(&B);
 92:   PetscFinalize();
 93:   return 0;
 94: }

 96: /*TEST

 98:     test:
 99:       nsize: {{1 2 3 4}}

101: TEST*/