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