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