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