Actual source code: ex229.c

  1: static char help[] = "Test MATMFFD for the rectangular case\n\n";

  3: #include <petscmat.h>

  5: static PetscErrorCode myF(void *ctx, Vec x, Vec y)
  6: {
  7:   const PetscScalar *ax;
  8:   PetscScalar       *ay;
  9:   PetscInt           i, j, m, n;

 11:   PetscFunctionBegin;
 12:   PetscCall(VecGetArrayRead(x, &ax));
 13:   PetscCall(VecGetArray(y, &ay));
 14:   PetscCall(VecGetLocalSize(y, &m));
 15:   PetscCall(VecGetLocalSize(x, &n));
 16:   for (i = 0; i < m; i++) {
 17:     PetscScalar xx, yy;

 19:     yy = 0.0;
 20:     for (j = 0; j < n; j++) {
 21:       xx = PetscPowScalarInt(ax[j], i + 1);
 22:       yy += xx;
 23:     }
 24:     ay[i] = yy;
 25:   }
 26:   PetscCall(VecRestoreArray(y, &ay));
 27:   PetscCall(VecRestoreArrayRead(x, &ax));
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

 31: int main(int argc, char **args)
 32: {
 33:   Mat      A, B;
 34:   Vec      base;
 35:   PetscInt m = 3, n = 2;

 37:   PetscFunctionBeginUser;
 38:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 39:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 40:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 41:   PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, n, &A));
 42:   PetscCall(MatCreateVecs(A, &base, NULL));
 43:   PetscCall(VecSet(base, 2.0));
 44:   PetscCall(MatMFFDSetFunction(A, myF, NULL));
 45:   PetscCall(MatMFFDSetBase(A, base, NULL));
 46:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 47:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 48:   PetscCall(MatComputeOperator(A, NULL, &B));
 49:   PetscCall(VecDestroy(&base));
 50:   PetscCall(MatDestroy(&A));
 51:   PetscCall(MatDestroy(&B));
 52:   PetscCall(PetscFinalize());
 53:   return 0;
 54: }

 56: /*TEST

 58:     test:
 59:       nsize: {{1 2 3 4}}

 61: TEST*/