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 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 38 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 39 PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); 40 PetscCall(MatCreateMFFD(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m,n,&A)); 41 PetscCall(MatCreateVecs(A,&base,NULL)); 42 PetscCall(VecSet(base,2.0)); 43 PetscCall(MatMFFDSetFunction(A,myF,NULL)); 44 PetscCall(MatMFFDSetBase(A,base,NULL)); 45 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 46 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 47 PetscCall(MatComputeOperator(A,NULL,&B)); 48 PetscCall(VecDestroy(&base)); 49 PetscCall(MatDestroy(&A)); 50 PetscCall(MatDestroy(&B)); 51 PetscCall(PetscFinalize()); 52 return 0; 53 } 54 55 /*TEST 56 57 test: 58 nsize: {{1 2 3 4}} 59 60 TEST*/ 61