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";
4
5 #include <petscmat.h>
6 #include <petscksp.h>
7
main(int argc,char ** args)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;
18
19 PetscFunctionBeginUser;
20 PetscCall(PetscInitialize(&argc, &args, NULL, help));
21
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));
31
32 /*
33 If the load matrix is larger than 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 }
51
52 /* Set up solution */
53 PetscCall(VecDuplicate(b, &x));
54 PetscCall(VecDuplicate(b, &u));
55 PetscCall(VecSet(x, 0.0));
56
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());
65
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));
73
74 /* Cleanup */
75 PetscCall(KSPDestroy(&ksp));
76 PetscCall(VecDestroy(&x));
77 PetscCall(VecDestroy(&b));
78 PetscCall(VecDestroy(&u));
79 PetscCall(MatDestroy(&A));
80
81 PetscCall(PetscFinalize());
82 return 0;
83 }
84
85 /*TEST
86
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)
90
91 TEST*/
92