2: static char help[] = "Tests the aSA multigrid code.\n"
  3: "Parameters:\n"
  4: "-n n          to use a matrix size of n\n";

  6: #include <petscdmda.h>
  7: #include <petscksp.h>
  8: #include <petscpcasa.h>

 10: PetscErrorCode  Create1dLaplacian(PetscInt,Mat*);
 11: PetscErrorCode  CalculateRhs(Vec);

 15: int main(int Argc,char **Args)
 16: {
 17:   PetscInt        n = 60;
 18:   PetscErrorCode  ierr;
 19:   Mat             cmat;
 20:   Vec             b,x;
 21:   KSP             kspmg;
 22:   PC              pcmg;
 23:   DM              da;

 25:   PetscInitialize(&Argc,&Args,(char *)0,help);
 26: 
 27:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 28:   Create1dLaplacian(n,&cmat);
 29:   MatGetVecs(cmat,&x,0);
 30:   MatGetVecs(cmat,&b,0);
 31:   CalculateRhs(b);
 32:   VecSet(x,0.0);

 34:   KSPCreate(PETSC_COMM_WORLD,&kspmg);
 35:   KSPSetType(kspmg, KSPCG);

 37:   KSPGetPC(kspmg,&pcmg);
 38:   PCSetType(pcmg,PCASA);

 40:   /* maybe user wants to override some of the choices */
 41:   KSPSetFromOptions(kspmg);

 43:   KSPSetOperators(kspmg,cmat,cmat,DIFFERENT_NONZERO_PATTERN);

 45:   DMDACreate1d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, n, 1, 1, 0, &da);
 46:   DMDASetRefinementFactor(da, 3, 3, 3);
 47:   PCSetDM(pcmg, (DM) da);

 49:   PCASASetTolerances(pcmg, 1.e-10, 1.e-10, PETSC_DEFAULT,PETSC_DEFAULT);

 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 52:                       Solve the linear system
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 55:   KSPSolve(kspmg,b,x);
 56:   KSPDestroy(&kspmg);
 57:   VecDestroy(&x);
 58:   VecDestroy(&b);
 59:   MatDestroy(&cmat);
 60:   DMDestroy(&da);
 61:   PetscFinalize();
 62:   return 0;
 63: }

 65: /* --------------------------------------------------------------------- */
 68: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
 69: {
 70:   PetscScalar    mone = -1.0,two = 2.0;
 71:   PetscInt       i,j,loc_start,loc_end;

 75:   MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n,PETSC_DECIDE, PETSC_NULL, PETSC_DECIDE, PETSC_NULL, mat);

 77:   MatGetOwnershipRange(*mat,&loc_start,&loc_end);
 78:   for (i=loc_start; i<loc_end; i++) {
 79:     if(i>0)   { j=i-1; MatSetValues(*mat,1,&i,1,&j,&mone,INSERT_VALUES); }
 80:     MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
 81:     if(i<n-1) { j=i+1; MatSetValues(*mat,1,&i,1,&j,&mone,INSERT_VALUES); }
 82:   }
 83:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
 84:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
 85:   return(0);
 86: }
 87: /* --------------------------------------------------------------------- */
 90: PetscErrorCode CalculateRhs(Vec u)
 91: {
 93:   PetscInt       i,n,loc_start,loc_end;
 94:   PetscReal      h;
 95:   PetscScalar    uu;

 98:   VecGetSize(u,&n);
 99:   VecGetOwnershipRange(u,&loc_start,&loc_end);
100:   h = 1.0/((PetscReal)(n+1));
101:   uu = 2.0*h*h;
102:   for (i=loc_start; i<loc_end; i++) {
103:     VecSetValues(u,1,&i,&uu,INSERT_VALUES);
104:   }
105:   VecAssemblyBegin(u);
106:   VecAssemblyEnd(u);

108:   return(0);
109: }