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: }