Actual source code: ex35.c

petsc-3.4.2 2013-07-02
  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 <petscdmda.h>
 12: #include <petscksp.h>

 14: extern PetscErrorCode ComputeMatrix(DM,Mat);
 15: extern PetscErrorCode ComputeRHS(DM,Vec);

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

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

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

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

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

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

 68:   KSPSolve(ksp,b,x);

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

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