Actual source code: ex18.c
1: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
2: Input arguments are:\n\
3: -f <input_file> : file to load. For example see $PETSC_DIR/share/petsc/datafiles/matrices\n\n";
5: #include <petscmat.h>
6: #include <petscksp.h>
8: int main(int argc, char **args)
9: {
10: PetscInt its, m, n, mvec;
11: PetscReal norm;
12: Vec x, b, u;
13: Mat A;
14: KSP ksp;
15: char file[PETSC_MAX_PATH_LEN];
16: PetscViewer fd;
17: PetscLogStage stage1;
19: PetscFunctionBeginUser;
20: PetscCall(PetscInitialize(&argc, &args, NULL, help));
22: /* Read matrix and RHS */
23: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), NULL));
24: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
25: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
26: PetscCall(MatSetType(A, MATSEQAIJ));
27: PetscCall(MatLoad(A, fd));
28: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
29: PetscCall(VecLoad(b, fd));
30: PetscCall(PetscViewerDestroy(&fd));
32: /*
33: If the load matrix is larger then the vector, due to being padded
34: to match the blocksize then create a new padded vector
35: */
36: PetscCall(MatGetSize(A, &m, &n));
37: PetscCall(VecGetSize(b, &mvec));
38: if (m > mvec) {
39: Vec tmp;
40: PetscScalar *bold, *bnew;
41: /* create a new vector b by padding the old one */
42: PetscCall(VecCreate(PETSC_COMM_WORLD, &tmp));
43: PetscCall(VecSetSizes(tmp, PETSC_DECIDE, m));
44: PetscCall(VecSetFromOptions(tmp));
45: PetscCall(VecGetArray(tmp, &bnew));
46: PetscCall(VecGetArray(b, &bold));
47: PetscCall(PetscArraycpy(bnew, bold, mvec));
48: PetscCall(VecDestroy(&b));
49: b = tmp;
50: }
52: /* Set up solution */
53: PetscCall(VecDuplicate(b, &x));
54: PetscCall(VecDuplicate(b, &u));
55: PetscCall(VecSet(x, 0.0));
57: /* Solve system */
58: PetscCall(PetscLogStageRegister("Stage 1", &stage1));
59: PetscCall(PetscLogStagePush(stage1));
60: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
61: PetscCall(KSPSetOperators(ksp, A, A));
62: PetscCall(KSPSetFromOptions(ksp));
63: PetscCall(KSPSolve(ksp, b, x));
64: PetscCall(PetscLogStagePop());
66: /* Show result */
67: PetscCall(MatMult(A, x, u));
68: PetscCall(VecAXPY(u, -1.0, b));
69: PetscCall(VecNorm(u, NORM_2, &norm));
70: PetscCall(KSPGetIterationNumber(ksp, &its));
71: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of iterations = %3" PetscInt_FMT "\n", its));
72: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm));
74: /* Cleanup */
75: PetscCall(KSPDestroy(&ksp));
76: PetscCall(VecDestroy(&x));
77: PetscCall(VecDestroy(&b));
78: PetscCall(VecDestroy(&u));
79: PetscCall(MatDestroy(&A));
81: PetscCall(PetscFinalize());
82: return 0;
83: }
85: /*TEST
87: test:
88: args: -ksp_gmres_cgs_refinement_type refine_always -f ${DATAFILESPATH}/matrices/arco1 -ksp_monitor_short
89: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
91: TEST*/