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 { 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 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, (char *)0, 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 61 TEST*/ 62