Actual source code: ex17.c
2: static char help[] = "Tests the use of MatSolveTranspose().\n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat C, A;
9: PetscInt i, j, m = 5, n = 5, Ii, J;
10: PetscScalar v, five = 5.0, one = 1.0;
11: IS isrow, row, col;
12: Vec x, u, b;
13: PetscReal norm;
14: MatFactorInfo info;
17: PetscInitialize(&argc, &args, (char *)0, help);
18: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
19: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
21: MatCreateSeqAIJ(PETSC_COMM_SELF, m * n, m * n, 5, NULL, &C);
22: MatSetUp(C);
24: /* create the matrix for the five point stencil, YET AGAIN*/
25: for (i = 0; i < m; i++) {
26: for (j = 0; j < n; j++) {
27: v = -1.0;
28: Ii = j + n * i;
29: if (i > 0) {
30: J = Ii - n;
31: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
32: }
33: if (i < m - 1) {
34: J = Ii + n;
35: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
36: }
37: if (j > 0) {
38: J = Ii - 1;
39: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
40: }
41: if (j < n - 1) {
42: J = Ii + 1;
43: MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES);
44: }
45: v = 4.0;
46: MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
47: }
48: }
49: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
50: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
52: ISCreateStride(PETSC_COMM_SELF, (m * n) / 2, 0, 2, &isrow);
53: MatZeroRowsIS(C, isrow, five, 0, 0);
55: VecCreateSeq(PETSC_COMM_SELF, m * n, &u);
56: VecDuplicate(u, &x);
57: VecDuplicate(u, &b);
58: VecSet(u, one);
60: MatMultTranspose(C, u, b);
62: /* Set default ordering to be Quotient Minimum Degree; also read
63: orderings from the options database */
64: MatGetOrdering(C, MATORDERINGQMD, &row, &col);
66: MatFactorInfoInitialize(&info);
67: MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A);
68: MatLUFactorSymbolic(A, C, row, col, &info);
69: MatLUFactorNumeric(A, C, &info);
70: MatSolveTranspose(A, b, x);
72: ISView(row, PETSC_VIEWER_STDOUT_SELF);
73: VecAXPY(x, -1.0, u);
74: VecNorm(x, NORM_2, &norm);
75: PetscPrintf(PETSC_COMM_SELF, "Norm of error %g\n", (double)norm);
77: ISDestroy(&row);
78: ISDestroy(&col);
79: ISDestroy(&isrow);
80: VecDestroy(&u);
81: VecDestroy(&x);
82: VecDestroy(&b);
83: MatDestroy(&C);
84: MatDestroy(&A);
85: PetscFinalize();
86: return 0;
87: }
89: /*TEST
91: test:
93: TEST*/