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:   VecGetArrayRead(x, &ax);
 12:   VecGetArray(y, &ay);
 13:   VecGetLocalSize(y, &m);
 14:   VecGetLocalSize(x, &n);
 15:   for (i = 0; i < m; i++) {
 16:     PetscScalar xx, yy;

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

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

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

 55: /*TEST

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

 60: TEST*/