Actual source code: ex27.c
1: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
2: Test MatMatSolve(). Input parameters include\n\
3: -f <input_file> : file to load \n\n";
5: /*
6: Usage:
7: ex27 -f0 <mat_binaryfile>
8: */
10: #include <petscksp.h>
11: extern PetscErrorCode PCShellApply_Matinv(PC, Vec, Vec);
13: int main(int argc, char **args)
14: {
15: KSP ksp;
16: Mat A, B, F, X;
17: Vec x, b, u; /* approx solution, RHS, exact solution */
18: PetscViewer fd; /* viewer */
19: char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
20: PetscBool flg;
21: PetscInt M, N, i, its;
22: PetscReal norm;
23: PetscScalar val = 1.0;
24: PetscMPIInt size;
25: PC pc;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &args, NULL, help));
29: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
30: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
32: /* Read matrix and right-hand-side vector */
33: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg));
34: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
36: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd));
37: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
38: PetscCall(MatSetType(A, MATAIJ));
39: PetscCall(MatLoad(A, fd));
40: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
41: PetscCall(VecLoad(b, fd));
42: PetscCall(PetscViewerDestroy(&fd));
44: /*
45: If the loaded matrix is larger than the vector (due to being padded
46: to match the block size of the system), then create a new padded vector.
47: */
48: {
49: PetscInt m, n, j, mvec, start, end, indx;
50: Vec tmp;
51: PetscScalar *bold;
53: /* Create a new vector b by padding the old one */
54: PetscCall(MatGetLocalSize(A, &m, &n));
55: PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);
56: PetscCall(VecCreate(PETSC_COMM_WORLD, &tmp));
57: PetscCall(VecSetSizes(tmp, m, PETSC_DECIDE));
58: PetscCall(VecSetFromOptions(tmp));
59: PetscCall(VecGetOwnershipRange(b, &start, &end));
60: PetscCall(VecGetLocalSize(b, &mvec));
61: PetscCall(VecGetArray(b, &bold));
62: for (j = 0; j < mvec; j++) {
63: indx = start + j;
64: PetscCall(VecSetValues(tmp, 1, &indx, bold + j, INSERT_VALUES));
65: }
66: PetscCall(VecRestoreArray(b, &bold));
67: PetscCall(VecDestroy(&b));
68: PetscCall(VecAssemblyBegin(tmp));
69: PetscCall(VecAssemblyEnd(tmp));
70: b = tmp;
71: }
72: PetscCall(VecDuplicate(b, &x));
73: PetscCall(VecDuplicate(b, &u));
74: PetscCall(VecSet(x, 0.0));
76: /* Create dense matrices B and X. Set B as an identity matrix */
77: PetscCall(MatGetSize(A, &M, &N));
78: PetscCall(MatCreate(MPI_COMM_SELF, &B));
79: PetscCall(MatSetSizes(B, M, N, M, N));
80: PetscCall(MatSetType(B, MATSEQDENSE));
81: PetscCall(MatSeqDenseSetPreallocation(B, NULL));
82: for (i = 0; i < M; i++) PetscCall(MatSetValues(B, 1, &i, 1, &i, &val, INSERT_VALUES));
83: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
84: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
86: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X));
88: /* Compute X=inv(A) by MatMatSolve() */
89: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
90: PetscCall(KSPSetOperators(ksp, A, A));
91: PetscCall(KSPGetPC(ksp, &pc));
92: PetscCall(PCSetType(pc, PCLU));
93: PetscCall(KSPSetFromOptions(ksp));
94: PetscCall(KSPSetUp(ksp));
95: PetscCall(PCFactorGetMatrix(pc, &F));
96: PetscCall(MatMatSolve(F, B, X));
97: PetscCall(MatDestroy(&B));
99: /* Now, set X=inv(A) as a preconditioner */
100: PetscCall(PCSetType(pc, PCSHELL));
101: PetscCall(PCShellSetContext(pc, X));
102: PetscCall(PCShellSetApply(pc, PCShellApply_Matinv));
103: PetscCall(KSPSetFromOptions(ksp));
105: /* Solve preconditioned system A*x = b */
106: PetscCall(KSPSolve(ksp, b, x));
107: PetscCall(KSPGetIterationNumber(ksp, &its));
109: /* Check error */
110: PetscCall(MatMult(A, x, u));
111: PetscCall(VecAXPY(u, -1.0, b));
112: PetscCall(VecNorm(u, NORM_2, &norm));
113: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));
116: /* Free work space. */
117: PetscCall(MatDestroy(&X));
118: PetscCall(MatDestroy(&A));
119: PetscCall(VecDestroy(&b));
120: PetscCall(VecDestroy(&u));
121: PetscCall(VecDestroy(&x));
122: PetscCall(KSPDestroy(&ksp));
123: PetscCall(PetscFinalize());
124: return 0;
125: }
127: PetscErrorCode PCShellApply_Matinv(PC pc, Vec xin, Vec xout)
128: {
129: Mat X;
131: PetscFunctionBeginUser;
132: PetscCall(PCShellGetContext(pc, &X));
133: PetscCall(MatMult(X, xin, xout));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: /*TEST
139: test:
140: args: -f ${DATAFILESPATH}/matrices/small
141: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
142: output_file: output/ex27.out
144: TEST*/