1 static char help[] = "Test MATMFFD for the rectangular case\n\n";
2
3 #include <petscmat.h>
4
myF(PetscCtx ctx,Vec x,Vec y)5 static PetscErrorCode myF(PetscCtx ctx, Vec x, Vec y)
6 {
7 const PetscScalar *ax;
8 PetscScalar *ay;
9 PetscInt i, j, m, n;
10
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;
18
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 }
30
main(int argc,char ** args)31 int main(int argc, char **args)
32 {
33 Mat A, B;
34 Vec base;
35 PetscInt m = 3, n = 2;
36
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 }
55
56 /*TEST
57
58 test:
59 nsize: {{1 2 3 4}}
60 output_file: output/empty.out
61
62 TEST*/
63