Actual source code: ex35.c

petsc-3.3-p6 2013-02-11
  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>
 13: #include <petscpcmg.h>

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

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

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

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

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

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

 63:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 64:   KSPSetFromOptions(ksp);
 65:   KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
 66:   KSPGetPC(ksp,&pc);
 67:   PCSetDM(pc,(DM)da);
 68: 
 69:   KSPSolve(ksp,b,x);

 71:   /* check final residual */
 72:   flg  = PETSC_FALSE;
 73:   PetscOptionsGetBool(PETSC_NULL, "-check_final_residual", &flg,PETSC_NULL);
 74:   if (flg){
 75:     Vec            b1;
 76:     PetscReal      norm;
 77:     KSPGetSolution(ksp,&x);
 78:     VecDuplicate(b,&b1);
 79:     MatMult(A,x,b1);
 80:     VecAXPY(b1,-1.0,b);
 81:     VecNorm(b1,NORM_2,&norm);
 82:     PetscPrintf(PETSC_COMM_WORLD,"Final residual %g\n",norm);
 83:     VecDestroy(&b1);
 84:   }
 85: 
 86:   KSPDestroy(&ksp);
 87:   VecDestroy(&x);
 88:   VecDestroy(&b);
 89:   MatDestroy(&A);
 90:   DMDestroy(&da);
 91:   PetscFinalize();
 92:   return 0;
 93: }