Actual source code: ex50.c


  2: static char help[] = "Tests point block Jacobi and ILU for different block sizes\n\n";

  4: #include <petscksp.h>

  6: int main(int argc, char **args)
  7: {
  8:   Vec         x, b, u;
  9:   Mat         A;    /* linear system matrix */
 10:   KSP         ksp;  /* linear solver context */
 11:   PetscRandom rctx; /* random number generator context */
 12:   PetscReal   norm; /* norm of solution error */
 13:   PetscInt    i, j, k, l, n = 27, its, bs = 2, Ii, J;
 14:   PetscScalar v;

 17:   PetscInitialize(&argc, &args, (char *)0, help);
 18:   PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL);
 19:   PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);

 21:   MatCreate(PETSC_COMM_WORLD, &A);
 22:   MatSetSizes(A, n * bs, n * bs, PETSC_DETERMINE, PETSC_DETERMINE);
 23:   MatSetBlockSize(A, bs);
 24:   MatSetFromOptions(A);
 25:   MatSetUp(A);

 27:   /*
 28:      Don't bother to preallocate matrix
 29:   */
 30:   PetscRandomCreate(PETSC_COMM_SELF, &rctx);
 31:   for (i = 0; i < n; i++) {
 32:     for (j = 0; j < n; j++) {
 33:       PetscRandomGetValue(rctx, &v);
 34:       if (PetscRealPart(v) < .25 || i == j) {
 35:         for (k = 0; k < bs; k++) {
 36:           for (l = 0; l < bs; l++) {
 37:             PetscRandomGetValue(rctx, &v);
 38:             Ii = i * bs + k;
 39:             J  = j * bs + l;
 40:             if (Ii == J) v += 10.;
 41:             MatSetValue(A, Ii, J, v, INSERT_VALUES);
 42:           }
 43:         }
 44:       }
 45:     }
 46:   }

 48:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 49:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 50:   MatCreateVecs(A, &u, &b);
 51:   VecDuplicate(u, &x);
 52:   VecSet(u, 1.0);
 53:   MatMult(A, u, b);

 55:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:                 Create the linear solver and set various options
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 59:   /*
 60:      Create linear solver context
 61:   */
 62:   KSPCreate(PETSC_COMM_WORLD, &ksp);

 64:   /*
 65:      Set operators. Here the matrix that defines the linear system
 66:      also serves as the preconditioning matrix.
 67:   */
 68:   KSPSetOperators(ksp, A, A);

 70:   KSPSetFromOptions(ksp);

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:                       Solve the linear system
 74:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 76:   KSPSolve(ksp, b, x);

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:                       Check solution and clean up
 80:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 82:   /*
 83:      Check the error
 84:   */
 85:   VecAXPY(x, -1.0, u);
 86:   VecNorm(x, NORM_2, &norm);
 87:   KSPGetIterationNumber(ksp, &its);

 89:   /*
 90:      Print convergence information.  PetscPrintf() produces a single
 91:      print statement from all processes that share a communicator.
 92:      An alternative is PetscFPrintf(), which prints to a file.
 93:   */
 94:   if (norm > .1) PetscPrintf(PETSC_COMM_WORLD, "Norm of residual %g iterations %" PetscInt_FMT " bs %" PetscInt_FMT "\n", (double)norm, its, bs);

 96:   /*
 97:      Free work space.  All PETSc objects should be destroyed when they
 98:      are no longer needed.
 99:   */
100:   KSPDestroy(&ksp);
101:   VecDestroy(&u);
102:   VecDestroy(&x);
103:   VecDestroy(&b);
104:   MatDestroy(&A);
105:   PetscRandomDestroy(&rctx);

107:   /*
108:      Always call PetscFinalize() before exiting a program.  This routine
109:        - finalizes the PETSc libraries as well as MPI
110:        - provides summary and diagnostic information if certain runtime
111:          options are chosen (e.g., -log_view).
112:   */
113:   PetscFinalize();
114:   return 0;
115: }

117: /*TEST
118:   testset:
119:     args: -mat_type {{aij baij}}

121:     test:
122:       args: -bs {{1 2 3 4 5 6 7 8}} -pc_type pbjacobi

124:     test:
125:       suffix: 2
126:       args: -bs {{8 9 10 11 12 13 14 15}} -pc_type ilu

128: TEST*/