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(0); 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