Actual source code: ex60.c

  1: static char help[] = "Working out corner cases of the ASM preconditioner.\n\n";

  3: #include <petscksp.h>

  5: int main(int argc,char **args)
  6: {
  7:   KSP            ksp;
  8:   PC             pc;
  9:   Mat            A;
 10:   Vec            u, x, b;
 11:   PetscReal      error;
 12:   PetscMPIInt    rank, size, sized;
 13:   PetscInt       M = 8, N = 8, m, n, rstart, rend, r;
 14:   PetscBool      userSubdomains = PETSC_FALSE;

 16:   PetscInitialize(&argc, &args, NULL,help);
 17:   PetscOptionsGetInt(NULL,NULL, "-M", &M, NULL);
 18:   PetscOptionsGetInt(NULL,NULL, "-N", &N, NULL);
 19:   PetscOptionsGetBool(NULL,NULL, "-user_subdomains", &userSubdomains, NULL);
 20:   /* Do parallel decomposition */
 21:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 22:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 23:   sized = (PetscMPIInt) PetscSqrtReal((PetscReal) size);
 27:   /* Assemble the matrix for the five point stencil, YET AGAIN
 28:        Every other process will be empty */
 29:   MatCreate(PETSC_COMM_WORLD, &A);
 30:   m    = (sized > 1) ? (rank % 2) ? 0 : 2*M/sized : M;
 31:   n    = N/sized;
 32:   MatSetSizes(A, m*n, m*n, M*N, M*N);
 33:   MatSetFromOptions(A);
 34:   MatSetUp(A);
 35:   MatGetOwnershipRange(A, &rstart, &rend);
 36:   for (r = rstart; r < rend; ++r) {
 37:     const PetscScalar diag = 4.0, offdiag = -1.0;
 38:     const PetscInt    i    = r/N;
 39:     const PetscInt    j    = r - i*N;
 40:     PetscInt          c;

 42:     if (i > 0)   {c = r - n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
 43:     if (i < M-1) {c = r + n; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
 44:     if (j > 0)   {c = r - 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
 45:     if (j < N-1) {c = r + 1; MatSetValues(A, 1, &r, 1, &c, &offdiag, INSERT_VALUES);}
 46:     MatSetValues(A, 1, &r, 1, &r, &diag, INSERT_VALUES);
 47:   }
 48:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 50:   /* Setup Solve */
 51:   MatCreateVecs(A, &x, &b);
 52:   VecDuplicate(x, &u);
 53:   VecSet(u, 1.0);
 54:   MatMult(A, u, b);
 55:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 56:   KSPSetOperators(ksp, A, A);
 57:   KSPGetPC(ksp, &pc);
 58:   PCSetType(pc, PCASM);
 59:   /* Setup ASM by hand */
 60:   if (userSubdomains) {
 61:     IS        is;
 62:     PetscInt *rows;

 64:     /* Use no overlap for now */
 65:     PetscMalloc1(rend-rstart, &rows);
 66:     for (r = rstart; r < rend; ++r) rows[r-rstart] = r;
 67:     ISCreateGeneral(PETSC_COMM_SELF, rend-rstart, rows, PETSC_OWN_POINTER, &is);
 68:     PCASMSetLocalSubdomains(pc, 1, &is, &is);
 69:     ISDestroy(&is);
 70:   }
 71:   KSPSetFromOptions(ksp);
 72:   /* Solve and Compare */
 73:   KSPSolve(ksp, b, x);
 74:   VecAXPY(x, -1.0, u);
 75:   VecNorm(x, NORM_INFINITY, &error);
 76:   PetscPrintf(PETSC_COMM_WORLD, "Infinity norm of the error: %g\n", (double) error);
 77:   /* Cleanup */
 78:   KSPDestroy(&ksp);
 79:   MatDestroy(&A);
 80:   VecDestroy(&u);
 81:   VecDestroy(&x);
 82:   VecDestroy(&b);
 83:   PetscFinalize();
 84:   return 0;
 85: }

 87: /*TEST

 89:    test:
 90:       suffix: 0
 91:       args: -ksp_view

 93:    test:
 94:       requires: !sycl kokkos_kernels
 95:       suffix: 0_kokkos
 96:       args: -ksp_view -mat_type aijkokkos

 98:    test:
 99:       requires: cuda
100:       suffix: 0_cuda
101:       args: -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

103:    test:
104:       suffix: 1
105:       nsize: 4
106:       args: -ksp_view

108:    test:
109:       requires: !sycl kokkos_kernels
110:       suffix: 1_kokkos
111:       nsize: 4
112:       args: -ksp_view -mat_type aijkokkos

114:    test:
115:       requires: cuda
116:       suffix: 1_cuda
117:       nsize: 4
118:       args: -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

120:    test:
121:       suffix: 2
122:       nsize: 4
123:       args: -user_subdomains -ksp_view

125:    test:
126:       requires: !sycl kokkos_kernels
127:       suffix: 2_kokkos
128:       nsize: 4
129:       args: -user_subdomains -ksp_view -mat_type aijkokkos

131:    test:
132:       requires: cuda
133:       suffix: 2_cuda
134:       nsize: 4
135:       args: -user_subdomains -ksp_view -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

137: TEST*/