1 static char help[] = "Test MATMFFD for the rectangular case\n\n"; 2 3 #include <petscmat.h> 4 5 static PetscErrorCode myF(void *ctx, Vec x, Vec y) { 6 const PetscScalar *ax; 7 PetscScalar *ay; 8 PetscInt i, j, m, n; 9 10 PetscFunctionBegin; 11 PetscCall(VecGetArrayRead(x, &ax)); 12 PetscCall(VecGetArray(y, &ay)); 13 PetscCall(VecGetLocalSize(y, &m)); 14 PetscCall(VecGetLocalSize(x, &n)); 15 for (i = 0; i < m; i++) { 16 PetscScalar xx, yy; 17 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 PetscCall(VecRestoreArray(y, &ay)); 26 PetscCall(VecRestoreArrayRead(x, &ax)); 27 PetscFunctionReturn(0); 28 } 29 30 int main(int argc, char **args) { 31 Mat A, B; 32 Vec base; 33 PetscInt m = 3, n = 2; 34 35 PetscFunctionBeginUser; 36 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 37 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 38 PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL)); 39 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, m, n, &A)); 40 PetscCall(MatCreateVecs(A, &base, NULL)); 41 PetscCall(VecSet(base, 2.0)); 42 PetscCall(MatMFFDSetFunction(A, myF, NULL)); 43 PetscCall(MatMFFDSetBase(A, base, NULL)); 44 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 45 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatComputeOperator(A, NULL, &B)); 47 PetscCall(VecDestroy(&base)); 48 PetscCall(MatDestroy(&A)); 49 PetscCall(MatDestroy(&B)); 50 PetscCall(PetscFinalize()); 51 return 0; 52 } 53 54 /*TEST 55 56 test: 57 nsize: {{1 2 3 4}} 58 59 TEST*/ 60