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;
29: PetscInitialize(&argc, &args, (char *)0, help);
30: MPI_Comm_size(PETSC_COMM_WORLD, &size);
33: /* Read matrix and right-hand-side vector */
34: PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg);
37: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd);
38: MatCreate(PETSC_COMM_WORLD, &A);
39: MatSetType(A, MATAIJ);
40: MatLoad(A, fd);
41: VecCreate(PETSC_COMM_WORLD, &b);
42: VecLoad(b, fd);
43: 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: MatGetLocalSize(A, &m, &n);
57: VecCreate(PETSC_COMM_WORLD, &tmp);
58: VecSetSizes(tmp, m, PETSC_DECIDE);
59: VecSetFromOptions(tmp);
60: VecGetOwnershipRange(b, &start, &end);
61: VecGetLocalSize(b, &mvec);
62: VecGetArray(b, &bold);
63: for (j = 0; j < mvec; j++) {
64: indx = start + j;
65: VecSetValues(tmp, 1, &indx, bold + j, INSERT_VALUES);
66: }
67: VecRestoreArray(b, &bold);
68: VecDestroy(&b);
69: VecAssemblyBegin(tmp);
70: VecAssemblyEnd(tmp);
71: b = tmp;
72: }
73: VecDuplicate(b, &x);
74: VecDuplicate(b, &u);
75: VecSet(x, 0.0);
77: /* Create dense matric B and X. Set B as an identity matrix */
78: MatGetSize(A, &M, &N);
79: MatCreate(MPI_COMM_SELF, &B);
80: MatSetSizes(B, M, N, M, N);
81: MatSetType(B, MATSEQDENSE);
82: MatSeqDenseSetPreallocation(B, NULL);
83: for (i = 0; i < M; i++) MatSetValues(B, 1, &i, 1, &i, &val, INSERT_VALUES);
84: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
87: MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &X);
89: /* Compute X=inv(A) by MatMatSolve() */
90: KSPCreate(PETSC_COMM_WORLD, &ksp);
91: KSPSetOperators(ksp, A, A);
92: KSPGetPC(ksp, &pc);
93: PCSetType(pc, PCLU);
94: KSPSetFromOptions(ksp);
95: KSPSetUp(ksp);
96: PCFactorGetMatrix(pc, &F);
97: MatMatSolve(F, B, X);
98: MatDestroy(&B);
100: /* Now, set X=inv(A) as a preconditioner */
101: PCSetType(pc, PCSHELL);
102: PCShellSetContext(pc, X);
103: PCShellSetApply(pc, PCShellApply_Matinv);
104: KSPSetFromOptions(ksp);
106: /* Solve preconditioned system A*x = b */
107: KSPSolve(ksp, b, x);
108: KSPGetIterationNumber(ksp, &its);
110: /* Check error */
111: MatMult(A, x, u);
112: VecAXPY(u, -1.0, b);
113: VecNorm(u, NORM_2, &norm);
114: PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its);
115: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm);
117: /* Free work space. */
118: MatDestroy(&X);
119: MatDestroy(&A);
120: VecDestroy(&b);
121: VecDestroy(&u);
122: VecDestroy(&x);
123: KSPDestroy(&ksp);
124: PetscFinalize();
125: return 0;
126: }
128: PetscErrorCode PCShellApply_Matinv(PC pc, Vec xin, Vec xout)
129: {
130: Mat X;
133: PCShellGetContext(pc, &X);
134: MatMult(X, xin, xout);
135: return 0;
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*/