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*/