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;
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
21: PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, m * n, m * n, 5, NULL, &C));
22: PetscCall(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: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
32: }
33: if (i < m - 1) {
34: J = Ii + n;
35: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
36: }
37: if (j > 0) {
38: J = Ii - 1;
39: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
40: }
41: if (j < n - 1) {
42: J = Ii + 1;
43: PetscCall(MatSetValues(C, 1, &Ii, 1, &J, &v, INSERT_VALUES));
44: }
45: v = 4.0;
46: PetscCall(MatSetValues(C, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
47: }
48: }
49: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
50: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
52: PetscCall(ISCreateStride(PETSC_COMM_SELF, (m * n) / 2, 0, 2, &isrow));
53: PetscCall(MatZeroRowsIS(C, isrow, five, 0, 0));
55: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m * n, &u));
56: PetscCall(VecDuplicate(u, &x));
57: PetscCall(VecDuplicate(u, &b));
58: PetscCall(VecSet(u, one));
60: PetscCall(MatMultTranspose(C, u, b));
62: /* Set default ordering to be Quotient Minimum Degree; also read
63: orderings from the options database */
64: PetscCall(MatGetOrdering(C, MATORDERINGQMD, &row, &col));
66: PetscCall(MatFactorInfoInitialize(&info));
67: PetscCall(MatGetFactor(C, MATSOLVERPETSC, MAT_FACTOR_LU, &A));
68: PetscCall(MatLUFactorSymbolic(A, C, row, col, &info));
69: PetscCall(MatLUFactorNumeric(A, C, &info));
70: PetscCall(MatSolveTranspose(A, b, x));
72: PetscCall(ISView(row, PETSC_VIEWER_STDOUT_SELF));
73: PetscCall(VecAXPY(x, -1.0, u));
74: PetscCall(VecNorm(x, NORM_2, &norm));
75: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of error %g\n", (double)norm));
77: PetscCall(ISDestroy(&row));
78: PetscCall(ISDestroy(&col));
79: PetscCall(ISDestroy(&isrow));
80: PetscCall(VecDestroy(&u));
81: PetscCall(VecDestroy(&x));
82: PetscCall(VecDestroy(&b));
83: PetscCall(MatDestroy(&C));
84: PetscCall(MatDestroy(&A));
85: PetscCall(PetscFinalize());
86: return 0;
87: }
89: /*TEST
91: test:
93: TEST*/