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