Actual source code: ex35.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: /*
  3:   Used for testing AIJ matrix with all zeros.
  4: */

  6: static char help[] = "Used for Solving a linear system where the matrix has all zeros.\n\n";
  7: /*
  8:  Example: mpiexec -n <np> ./ex35 -dof 2 -mat_view -check_final_residual
  9:  */

 11: #include <petscdm.h>
 12: #include <petscdmda.h>
 13: #include <petscksp.h>

 17: int main(int argc,char **argv)
 18: {
 20:   KSP            ksp;
 21:   PC             pc;
 22:   Vec            x,b;
 23:   DM             da;
 24:   Mat            A;
 25:   PetscInt       dof=1;
 26:   PetscBool      flg;
 27:   PetscScalar    zero=0.0;

 29:   PetscInitialize(&argc,&argv,(char*)0,help);
 30:   PetscOptionsGetInt(NULL,"-dof",&dof,NULL);

 32:   DMDACreate(PETSC_COMM_WORLD,&da);
 33:   DMDASetDim(da,3);
 34:   DMDASetBoundaryType(da,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE);
 35:   DMDASetStencilType(da,DMDA_STENCIL_STAR);
 36:   DMDASetSizes(da,3,3,3);
 37:   DMDASetNumProcs(da,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 38:   DMDASetDof(da,dof);
 39:   DMDASetStencilWidth(da,1);
 40:   DMDASetOwnershipRanges(da,NULL,NULL,NULL);
 41:   DMSetFromOptions(da);
 42:   DMSetUp(da);

 44:   DMCreateGlobalVector(da,&x);
 45:   DMCreateGlobalVector(da,&b);
 46:   DMSetMatType(da,MATAIJ);
 47:   DMCreateMatrix(da,&A);
 48:   VecSet(b,zero);

 50:   /* Test sbaij matrix */
 51:   flg  = PETSC_FALSE;
 52:   PetscOptionsGetBool(NULL,"-test_sbaij",&flg,NULL);
 53:   if (flg) {
 54:     Mat sA;
 55:     MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
 56:     MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA);
 57:     MatDestroy(&A);
 58:     A    = sA;
 59:   }

 61:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 62:   KSPSetFromOptions(ksp);
 63:   KSPSetOperators(ksp,A,A);
 64:   KSPGetPC(ksp,&pc);
 65:   PCSetDM(pc,(DM)da);

 67:   KSPSolve(ksp,b,x);

 69:   /* check final residual */
 70:   flg  = PETSC_FALSE;
 71:   PetscOptionsGetBool(NULL, "-check_final_residual", &flg,NULL);
 72:   if (flg) {
 73:     Vec       b1;
 74:     PetscReal norm;
 75:     KSPGetSolution(ksp,&x);
 76:     VecDuplicate(b,&b1);
 77:     MatMult(A,x,b1);
 78:     VecAXPY(b1,-1.0,b);
 79:     VecNorm(b1,NORM_2,&norm);
 80:     PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
 81:     VecDestroy(&b1);
 82:   }

 84:   KSPDestroy(&ksp);
 85:   VecDestroy(&x);
 86:   VecDestroy(&b);
 87:   MatDestroy(&A);
 88:   DMDestroy(&da);
 89:   PetscFinalize();
 90:   return 0;
 91: }