Actual source code: ex3.c

  1: static char help[] = "Tests relaxation for dense matrices.\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat         C;
  8:   Vec         u, x, b, e;
  9:   PetscInt    i, n = 10, midx[3];
 10:   PetscScalar v[3];
 11:   PetscReal   omega = 1.0, norm;

 13:   PetscFunctionBeginUser;
 14:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 15:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-omega", &omega, NULL));
 16:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 18:   PetscCall(MatCreate(PETSC_COMM_SELF, &C));
 19:   PetscCall(MatSetSizes(C, n, n, n, n));
 20:   PetscCall(MatSetType(C, MATSEQDENSE));
 21:   PetscCall(MatSetUp(C));
 22:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &b));
 23:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &x));
 24:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &u));
 25:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &e));
 26:   PetscCall(VecSet(u, 1.0));
 27:   PetscCall(VecSet(x, 0.0));

 29:   v[0] = -1.;
 30:   v[1] = 2.;
 31:   v[2] = -1.;
 32:   for (i = 1; i < n - 1; i++) {
 33:     midx[0] = i - 1;
 34:     midx[1] = i;
 35:     midx[2] = i + 1;
 36:     PetscCall(MatSetValues(C, 1, &i, 3, midx, v, INSERT_VALUES));
 37:   }
 38:   i       = 0;
 39:   midx[0] = 0;
 40:   midx[1] = 1;
 41:   v[0]    = 2.0;
 42:   v[1]    = -1.;
 43:   PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES));
 44:   i       = n - 1;
 45:   midx[0] = n - 2;
 46:   midx[1] = n - 1;
 47:   v[0]    = -1.0;
 48:   v[1]    = 2.;
 49:   PetscCall(MatSetValues(C, 1, &i, 2, midx, v, INSERT_VALUES));

 51:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
 52:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));

 54:   PetscCall(MatMult(C, u, b));

 56:   for (i = 0; i < n; i++) {
 57:     PetscCall(MatSOR(C, b, omega, SOR_FORWARD_SWEEP, 0.0, 1, 1, x));
 58:     PetscCall(VecWAXPY(e, -1.0, x, u));
 59:     PetscCall(VecNorm(e, NORM_2, &norm));
 60:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "2-norm of error %g\n", (double)norm));
 61:   }
 62:   PetscCall(MatDestroy(&C));
 63:   PetscCall(VecDestroy(&x));
 64:   PetscCall(VecDestroy(&b));
 65:   PetscCall(VecDestroy(&u));
 66:   PetscCall(VecDestroy(&e));
 67:   PetscCall(PetscFinalize());
 68:   return 0;
 69: }

 71: /*TEST

 73:    test:

 75: TEST*/